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The  critical-direction  theory  is  extended  to  include  nonlinear  elements  through  the 
introduction  of  the  describing  function  and  a general  definition  of  the  critical  perturbation 
radius.  The  approach  allows  the  calculation  of  the  nonlinear  Nyquist  robust  stability 
margin  for  systems  with  both  structured  and  unstructured  uncertainties,  hence  providing  a 
conclusive  test  for  robust  stability.  Finally,  analytic  expressions  for  finding  the 
minimizing  amplitude  are  obtained  for  circular  uncertain  value  sets. 

The  robustness  of  predictive  controllers  is  studied  for  systems  with  real 
parametric  uncertainties  that  belong  to  an  ellipsoidal  or  hyperbox  domain.  A parametric 
stability  margin  is  proposed  as  a quantitative  measure  of  robust  stability  as  well  as  a 
robust  controller  design  methodology  that  involves  evaluating  the  parametric  stability 
margin  over  a range  of  values  for  the  predictive  control  tuning  parameters.  In  addition, 
the  standard  predictive  control  algorithm  is  also  modified  to  include  a disturbance  model 
to  help  improve  regulatory  performance.  The  resulting  disturbance  predictive  control 


Vll 


demonstrated  improved  performance  on  linear  systems.  The  performance  on  nonlinear 
cases  is  not  satisfactory  suggesting  that  further  refinements,  such  as  adaptive  control 
schemes,  may  be  necessary  to  handle  such  cases. 

Two  support  tools  for  operating  e mulsion  p olymerization  reactors  are  designed 
and  implemented,  namely  a continuous  sampling  and  dilution  system  and  an  online 
conversion  estimator  based  on  a densitometer.  The  dilution  system  is  fully  automated  for 
data  acquisition  and  real-time  control.  A dynamic  first-principles  model  is  developed  and 
found  to  agree  with  experimental  data.  A new  approach  to  estimate  conversion  on-line  is 
developed  through  the  use  of  a densitometer  and  the  calculations  are  improved  by  the  use 
of  an  excess-volume  model.  Both  a one  and  two-parameter  excess-volume  models  are 
analyzed,  and  examples  are  presented  to  illustrate  the  effectiveness  of  these  models. 
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CHAPTER  1 
INTRODUCTION 

1.1  Motivation 

Uncertainty  is  present  in  any  physical  system.  This  uncertainty  directly  translates 
to  the  model  of  the  system  used  for  control  design,  and  is  most  often  in  the  form  of 
variations  in  model  parameters.  Control  systems  are  required  to  be  robust  (i.e.,  it 
functions  satisfactorily  under  these  uncertainties),  and  the  design  of  these  control 
schemes  is  known  as  robust  control.  One  important  aspect  of  the  robust  control  problem 
is  determining  if  a control  system  is  robustly  stable  with  respect  to  a given  set  of 
uncertainties. 

Robust  stability  has  been  studied  since  the  earliest  days  of  feedback  control.  The 
earliest  robust  stability  studies  focused  on  frequency  domain  methods  such  as  those  based 
on  Bode  plots  and  Nyquist  plots  (Nyquist,  1932)  which  resulted  in  the  gain  and  phase 
stability  margins  that  are  still  used  today.  With  the  advent  of  the  space  race  of  the 
1960’s,  the  focus  of  control  engineers  was  shifted  away  from  frequency  domain  robust 
stability  methods  to  the  field  of  optimal  control.  The  linear  quadratic  regulator  (LQG) 
design,  developed  during  this  era,  appeared  to  give  controllers  with  good  stability 
properties.  However  in  the  late  1970’s  it  was  found  that  LQGs  lost  their  stability 
guarantees  when  uncertainty  was  present  in  the  control  system. 

Due  to  this  robust  instability,  Hx  optimal  control  (Doyle,  1983;  Doyle  et  ai, 
1989;  Francis  and  Zames,  1984;  Safonov  and  Verma,  1985;  Zames,  1981)  was  introduced 
to  effectively  deal  with  the  robust  stability  and  control  issues.  The  theory  considers 
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unstructured  uncertainties  where  the  only  information  known  about  the  uncertainty  is  a 
norm  bound.  Normally  more  information  is  known  about  the  uncertainty  than  a simple 
norm  bound  and  other  robust  stability  analysis  methods  have  been  developed  that 
consider  these  structured  uncertainties. 

One  robust  stability  analysis  method  for  structure  uncertainties  is  the  critical- 
direction  theory  developed  by  Latchman  and  Crisalle  (1995)  and  Baab  et  cil.  (2001) 
which  addresses  the  problem  of  robust  stability  of  systems  affected  by  uncertainties  that 
are  characterized  in  terms  of  arbitrary  frequency-domain  value  sets  that  can  be  convex  or 
non-convex.  The  critical  direction  theory  proposes  the  Nyquist  robust  stability  margin  as 
a measure  of  robust  stability  which  is  similar  to  other  techniques  that  have  been 
developed.  The  main  advantage  of  this  margin  is  that  determining  the  Nyquist  robust 
stability  margin  is  a tractable  problem  as  opposed  to  the  challenging  computations 
required  previously. 

Another  type  of  structured  uncertainty  is  real  parametric  uncertainty  in  the 
process  model  which  began  to  receive  renewed  attention  with  the  result  of  Kharitonov 
(1979)  on  the  stability  of  interval  polynomials.  This  result  allows  the  stability  analysis  of 
a linear  time  invariant  system  to  be  performed  as  real  uncertainty  is  added  to  the 
parameters  (Bhattacharyya  et  al.,  1995).  Consequently,  the  parametric  stability  margin  is 
defined  as  the  length  of  the  smallest  perturbation  in  the  parameters  which  just  destabilizes 
the  closed  loop  system.  The  parametric  stability  margin  is  useful  in  controller  design  as  a 
means  of  designing  robustly  stabilizing  controllers. 

Continuous  emulsion  polymerization  reactors  have  several  advantages  over  batch 
reactors;  among  them  are  (1)  high  reaction  rates  proving  high  yield,  high  molecular 
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weight  polymers,  (2)  easier  operation  and  control,  and  (3)  more  consistent  polymer 
quality.  However  there  are  major  problems  associated  with  the  operation  of  continuous 
reactors  such  as  no  true  steady  state,  agglomeration,  and  fouling.  The  control  of 
emulsion  polymerization  reactors  requires  advanced  control  strategies  like  those 
discussed  earlier  and  would  help  remove  or  reduce  these  problems.  However, 
subsystems  must  be  developed  that  allow  for  the  measurement  of  critical  properties  such 
as  novel  sensors  to  measure  conversion  and  particle  size  distribution. 

1.2  Objecive  and  Structure  of  Dissertation 

This  dissertation  is  divided  into  two  parts:  theoretical  process  control  and 
emulsion  polymerization  subsystems.  Theoretical  process  control  consists  of  three 
chapters  including  the  nonlinear  Nyquist  robust  stability  margin,  systematic  algorithms 
for  the  design  of  robust  predictive  controllers,  and  the  derivation  of  a disturbance 
predictive  controller  used  in  wastewater  treatment  plants.  The  emulsion  polymerization 
subsystems  include  the  determination  of  conversion  on-line  as  well  as  the  automation  of  a 
continuous  sampling  and  dilution  system  that  aids  spectroscopic  analysis. 

The  first  goal  of  this  dissertation  is  the  extension  of  the  critical  direction  theory  to 
the  more  general  case  of  nonlinear  systems  and  is  presented  in  Chapter  2.  The  key  to 
extending  the  theory  is  the  introduction  of  a generalized  definition  of  the  critical  point  in 
a fashion  that  preserves  all  previous  results.  Several  nonlinear  systems  are  analyzed  and 
the  minimizing  amplitude  discovered  that  maximizes  the  nonlinear  Nyquist  robust 
stability  margin.  The  generalized  critical  direction  theory  is  applied  to  nonlinear 
uncertain  systems,  and  is  used  to  calculate  the  required  nonlinear  Nyquist  robust  stability 
margin  with  high  precision  and  in  the  context  of  a computationally  manageable 
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Chapter  3 proposes  a systematic  algorithm  to  design  robust  predictive  controllers 
using  the  parametric  stability  margin.  Hyperbox  and  ellipsoidal  parametric  uncertainties 
are  considered  and  analyzed.  Finally,  examples  are  used  to  illustrate  the  design  of  the 
robust  predictive  controllers  for  both  types  of  uncertainty. 

The  last  section  of  theoretical  process  control  is  the  development  of  a disturbance 
predictive  controller  for  use  in  wastewater  treatment  plants  and  is  discussed  in  Chapter  4. 
An  brief  introduction  to  the  wastewater  treatment  plant  model  is  given  followed  by  the 
derivation  of  the  disturbance  predictive  controller.  The  controller  is  connected  to  both  a 
linear  and  nonlinear  model  of  the  wastewater  treatment  plant  and  compared  to  the 
predictive  controller  developed  by  Crisalle  et  al.  (1989). 

The  final  two  chapters  deal  with  subsystems  required  in  the  control  of  continuous 
emulsion  polymerization  reactors.  Chapter  5 introduces  an  excess  volume  model  to 
improve  conversion  estimates  from  densitometer  data.  The  numerical  issues  associated 
with  this  model  are  discussed  and  three  polymerization  examples  shown.  Chapter  6 
discusses  the  design  and  implementation  of  a continuous  sampling  and  dilution  system. 
The  schematic  and  control  scheme  for  the  dilution  system  is  shown  and  a dynamic  model 
derived.  Finally,  the  experimental  and  simulated  closed  loop  systems  are  compared  and 
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CHAPTER  2 

COMPUTATION  OF  THE  NONLINEAR  NYQUIST  ROBUST  STABILITY  MARGIN 

FOR  UNCERTAIN  SYSTEMS 

2.1.  Introduction 

Much  attention  has  been  given  to  the  problem  of  computing  a robust  stability 
margin  for  systems  with  frequency-domain  uncertainty  descriptions.  The  most  notable 
contributions  include  the  structured  singular  value  //  (Doyle,  1982)  and  the  multivariable 
stability  margin  km  (Safonov,  1982).  The  critical  direction  theory  developed  by 
Latchman  and  Crisalle  (1995)  and  by  Latchman  et  al.  (1997)  introduced  the  Nyquist 
robust  stability  margin  kN,  a scalar  measure  of  robustness  that  has  been  effectively 
computed  for  single-input/single-output  (SISO)  systems  with  both  convex  (Latchman  et 
al.,  1997)  and  non-convex  (Baab  et  al.,  2001)  uncertain  value  sets.  However,  all  previous 
work  has  dealt  with  linear  transfer  functions.  This  chapter  extends  the  critical  direction 
theory  to  nonlinear  systems,  while  preserving  the  key  definitions  and  notation  proposed 
for  the  linear  theory. 

2.2.  Generalization  of  the  Critical  Direction  Theory 
2.2.1  Preliminaries 

Consider  the  SISO  linear  time-invariant  system 

gO)  = go00  + £(s)  (2.1) 

where  g0(^)  is  a known  nominal  transfer  function  and  c>(s)eA  is  an  unknown 
perturbation  belonging  to  the  unstructured  uncertainty  family  A . Unstructured 
uncertainties  have  frequency-domain  value  sets  that  are  circles  with  known  radius  r(a> ) , 
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namely,  sets  of  the  form  A:=  j£(s):|<5(/<y)|<r(<y):r(<y)>0}.  The  uncertain  system 
(2.1)  is  arranged  in  a negative  feedback  configuration  featuring  a nonlinear  element  / ( e ) , 
as  shown  in  Figure  2.1.  For  simplicity  in  exposition,  in  Figure  2.1  we  introduce  a slight 
abuse  in  notation  given  that  the  map  f(e) : 91  -»  91  is  a time-domain  operator,  while 
g(s ) is  a Laplace-domain  operator. 


Figure  2.1.  Closed-loop  configuration  of  a nonlinear  operator  and  an 
uncertain  linear  plant  under  negative  feedback. 


The  robust  stability  analysis  proposed  involves  three  assumptions:  (i)  the  nominal 
plant  is  stable  under  unity  negative  feedback,  (ii)  the  nominal  system  g0(s)  and  the 
uncertain  system  g(s)  have  the  same  number  of  unstable  poles,  and  (iii)  the  higher 
harmonics  of  the  nonlinear  element  / ( e ) can  be  neglected. 

2.2.2  Brief  Review  of  the  Critical  Direction  Theory 

The  critical  direction  theory  was  initially  developed  for  linear  systems,  namely, 
for  the  special  case  of  Figure  2.1  where  f(e)  = 1 . Figure  2.2  shows  a typical  Nyquist 
diagram  illustrating  the  nominal  frequency  response  g0(ja>) . The  figure  also  shows  the 
point  representing  the  nominal  frequency  response  g0{ja> ) at  a specific  frequency  &>,, 
along  with  its  associated  circular  uncertainty  value  set  u(cox)  with  radius  r(a>).  The 
critical  perturbation  radius  pc((t)x ) for  the  circular  value-set  is  defined  as  the  distance 
from  the  nominal  point  g0(ja>\)  to  the  boundary  of  the  uncertain  circular  value  set,  i.e., 
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Pc(®)=H®)  i 


(2.2) 


Im 


Figure  2.2.  Critical  direction  theory  illustrating  a nominal  point  g0(jco ) 
and  its  associated  circular  value  set,  u(<yj) , and  critical  perturbation  radius 
pc(ry,). 

The  main  result  of  the  critical  direction  theory  for  linear  systems  is  summarized  in 
Theorem  2.1: 

Let  gn(jct >)  be  a nominal  SISO  system  with  an  associated  circular  uncertainty  set 
A . Then  the  closed  loop  system  given  in  Figure  1 with  f (e)  = 1 is  robustly  stable  with 
respect  to  uncertainties  S(s)  e A if  and  only  if 


The  proof  is  omitted  for  brevity.  Details  are  given  in  Latchman  (1997).  The  Nyquist 
robust-stability  margin  is  then  defined  for  SISO  systems  as 


l + g„0'«)| 


(2.3) 


Let  kN:=  max  kN(co) . Then  from  Theorem  2.1  it  is  obvious  that  a necessary  and 
sufficient  condition  for  robust  stability  is  given  by  the  inequality 
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kN  <\ 


(2.4) 


The  calculation  of  the  critical  perturbation  radius  is  often  the  most  c hallenging 


for  the  Nyquist  robust  stability  margin  (2.3);  however,  for  the  case  of  circular  value  sets, 
the  task  is  extremely  simple  given  the  relationship  (2.2). 

2.2.3  Brief  Review  of  the  Describing  Function  Method 

Consider  now  the  nonlinear  feedback  loop  shown  in  Figure  2.1  with /(e)  ^ 1. 

The  analysis  of  closed  loop  stability  hinges  on  determining  whether  the  limiting-stability 
condition  of  sustained  signal  oscillation  can  occur.  In  that  limiting  case,  the  input  e(t)  to 
the  nonlinear  function  /(e)  is  sinusoidal,  and  it  can  be  written  as  e{t)  = asmcot  where 

a>  0 is  the  signal  amplitude.  Using  a Fourier  series  expansion,  the  output  from  the 
nonlinear  function  can  be  written  in  the  form 


which  features  a mean  level  a0 , a set  of  fundamental  components  of  amplitudes  a\  and 
b] , and  higher  harmonic  components  of  amplitudes  ak  and  bk  at  frequencies 
km,  k = 2, 3,-", oo  (Khalil,  1992). 

The  describing  function  method  assumes  that  all  the  higher  harmonics  can  be 
neglected,  and  as  a consequence  the  nonlinear  function  / (e)  can  then  be  approximated 
by  its  describing  function  n{a) , also  referred  to  as  the  equivalent  gain 


problem  to  address  when  applying  the  critical  direction  theory  to  extract  numerical  values 


oo 


/ (a  sin  mt)  = a0  + ^ (ak  cos  kmt  + bk  sin  kmt) 


(2.5) 


a 


where 
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«i(®)  = - 


and 


The  describing  function  is  a complex  number  that  depends  on  the  amplitude  of  the 
sinusoid,  and  it  can  be  displayed  as  a locus  on  the  complex  plane. 

The  validity  o f this  method  relies  on  the  a ssumption  that  the  higher  harmonics 
produced  by  the  nonlinear  function  f(e ) can  be  discarded  (Leigh,  1983).  This  is  ensured 
in  all  cases  where  the  transfer  function  g(s)  behaves  as  a low-pass  filter. 


In  this  section  the  SISO  critical  direction  theory  is  extended  to  nonlinear  systems. 
The  concepts  of  critical  direction,  critical  perturbation  radius,  and  critical  point  are 
appropriately  redefined  to  take  into  account  the  nonlinearity  of  the  system. 

After  approximating  the  nonlinear  element  f (e)  by  its  describing  function  n(a) , 
the  characteristic  equation  for  the  control  configuration  in  Figure  2.1  is  of  the  form 


In  the  suite,  we  refer  to  the  Nyquist  map  of  -1/  n(a ) for  all  a > 0 as  the  critical 
loci.  The  robust  stability  of  the  closed  loop  is  ensured  if  the  value  sets  of  g(jco)  exclude 


2.3.  Critical  Direction  Theory  for  Nonlinear  Systems 


l + «(a)g(j)  = 0 


Rearranging  terms 


(2.6) 


all  points  belonging  to  the  critical  loci.  Under  the  assumptions  established  earlier,  the 
critical  loci  exclusion  principle  is  in  fact  a necessary  and  sufficient  condition  for  robust 
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closed  loop  stability.  The  following  definitions  introduce  concepts  illustrated  in  Figure 
2.3  and  Figure  2.4  that  are  useful  for  formulating  a nonlinear  critical  direction  theory  for 
systems  with  a convex  and  non-convex  uncertainty  description  respectively. 

1 .  The  critical  directions 


are  interpreted  as  a set  of  unit  vectors  originating  at  g0(jco)  and  pointing  towards  the 
critical  loci. 

2.  T he  c ritical  c one  9c(co)  maps  all  critical  directions  from  the  nominal  point  to  the 
critical  loci. 

3.  The  critical  template  uc(a> ) contains  the  uncertainty  defined  by  the  intersection  of  the 
critical  cone  and  uncertainty  value  set. 

4.  The  critical  perturbation  radii  for  convex  systems  is  defined  as 

pc(co,a)~  vcwt\a\z  = g0(j(D)  + adc(a>,a)  euc(ry)j  (2.7) 

«e9t+  3 

Equation  (2.7)  states  that  the  critical  perturbation  radius  for  a convex  critical  value  set 
uc(a))  is  the  distance  from  the  nominal  point  g0(co)  to  the  uncertain  value  set  boundary 
do(a>)  along  the  critical  direction.  Equation  (2.7),  however,  is  not  suitable  for  non- 
convex  critical  value  sets  oc(co) . A general  definition  of  the  critical  perturbation  radius 

similar  to  Baab  et  al.  (2001)  is  proposed  which  is  applicable  to  both  the  convex  and  non- 
convex  cases 


dc(jco,a)  = - 


g°U®)  + 1 + go{j(o)n{a)  \n(a)\ 


goU®)  + -TT 
n(a) 


1 | l + n(a) 


pc(co,a)\=- 


(2.8) 


n(a ) 


+ g0  ( jco ) +<^(a>)  otherwise 
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where 


%(co)  = min 
ze/3c((o) 


n{a ) 


■ + z 


(2.9) 


Im 


Figure  2.3.  A convex  value  set  with  critical  cone,  critical  perturbation 
radius,  and  critical  value  set  illustrated). 


Figure  2.4.  A non-convex  value  set  with  critical  cone,  critical  perturbation 
radius,  and  critical  value  set  illustrated. 
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Equation  (2.8)  represents  the  distance  from  the  critical  loci  to  the  point  in  /3c{co) 
that  is  closest  to  the  critical  loci.  The  upper  statement  in  equation  (2.8)  states  that  when 
the  critical  loci  is  not  an  element  of  the  uncertainty  value  set  v(co) , the  critical 
perturbation  radius  pc(a> ) is  defined  as  the  difference  between  the  two  distances, 
namely,  the  distance  from  the  critical  loci  to  the  nominal  point  (represented  by 


n{a) 


+ go{j(0) 


) and  the  distance  from  t he  critical  loci  to  the  closest  critical-boundary 


intersection  (represented  by  £(®)).  Conversely,  when  the  critical  loci  is  an  element  of 
the  uncertainty  value  set,  the  lower  statement  in  equation  (2.8)  states  that  the  critical 
perturbation  is  the  sum  of  the  two  respective  distances.  When  the  critical  value  set  is 
convex,  equation  (2.8)  reduces  to  equation  (2.7)  because  there  is  only  one  critical 
boundary  intersection.  To  calculate  the  critical  perturbation  radius  from  (2.8),  it  is 
necessary  to  have  full  knowledge  of  the  critical  boundary  intersections  /3c{a>)  and  to 
evaluate  whether  the  set  membership  condition  -1  /n(a)eu(co)  holds.  For  either 
definition,  it  should  be  noted  that  pc(co)  > 0 for  all  frequencies  and  is  a function  of 
frequency  and  of  the  amplitude  to  the  describing  function. 

At  a given  frequency®, , the  critical-line  set  is  defined  as  the  set  of  oriented  lines 
with  origin  at  the  nominal  point  g0{jcox),  and  that  pass  through  any  point  in  the  critical 
loci  -1  / n(a)  . T he  critical  perturbation  radius pc(cox, a)  i s i nterpreted  a s t he  d istance 
from  the  nominal  point  to  the  point  where  the  critical  line  defined  for  an  amplitude  a 
intersects  with  the  uncertain  value  set.  Finally,  the  Nyquist  robust  stability  margin  for 
nonlinear  systems  is  defined  as 
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kN(a>)\  = max 


pc(o),a) 


a>  0 


1 


n(a ) 


+ g0Ua) 


(2.10) 


or  equivalently, 


\n{d)\pc{a,a) 

kN(a>)  = max  | — — — — r 

a> o \g0(jco)  + n(a)\ 


(2.11) 


The  following  robust-stability  theorems  for  nonlinear  systems  use  the  same 
assumptions  invoked  previously  with  Theorem  2.2  illustrated  below. 

Let  ga(s)  be  a nominal  system  with  associated  circular  uncertainty  set  A . Then 
the  closed  loop  system  given  in  Figure  1 is  robustly  stable  with  respect  to  all 
uncertainties  S(s ) e A if  and  only  if 


kN(a> ) < 1 


(2.12) 


for  all  frequencies  co>  0 . 

Proof  A complete  proof  is  given  in  Latchman  and  Crisalle  (1997)  for  the  case  where 
uc(<y)  is  convex.  For  non-convex  systems,  the  general  definition  (2.8)  of  pc(co)  is 
exploited.  From  the  critical  loci  exclusion  principle  the  uncertain  closed  loop  is  stable  if 
and  only  if  -1  / n(a)  <£  u(co)  V co  . Therefore,  to  prove  that  (2.12)  is  sufficient  for  robust 
stability,  it  must  be  shown  that  if  kN  (co)  <1  V<y  then  -1  / n{a)  £ v(co)  . To  prove 
this  by  contradiction,  first  assume  that  kN(co)<l  \/a>  and  that  3 co  such  that 
-l  / n(a)  e v(a>) . Apply  definitions  (2.8)  and  (2.10)  for  a frequency  at  which 


-1  /n(a)ev(co)  to  yield 
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kN(co)  = 


pA°>) 

, +go(J<°) 

n(a) 

. , +g0Ua>) 
n(a) 

. +£o(» 

n(a) 

= 1 + 


(?(*) 


1 


n(a ) 


+ g0(J°)) 


where  £( co ) is  a non-negative  real  scalar  given  by  (2.9).  Therefore  &jy(<y)>l  for  at  least 
one  frequency  which  contradicts  the  assumption.  Thus,  if  kN(co)  < 1 Vco  then  it  must 
follow  that  -1  / n(a)  <£  o(a>)  \/co  . To  prove  that  (2.12)  is  necessary  for  robust  stability, 
it  must  be  shown  that  if  -l/n(a)  t v(co)  \/ a>  then  kN{co)  < 1 V co  . Apply  definitions 
(2.8)  and  (2.10)  if  -1  / «(a)  g v{co)  to  yield 


kN{co)  = 


pA®) 


i 


n(a ) 


+ g0Ua>) 


, +g0Ua) 

n(a) 

-A®) 

. +goU<») 

n(a) 

= 1- 


l 


n(a) 


+ go(Ja) 


V co 


where  %(a>)  is  given  by  (2.9).  Because  we  assumed  that  -1  / n(a)  £ o(a>)  it  follows  that 
-1  / n(a)  <£(3c(co)  and  thus  %{a>)  must  be  a positive  number.  This  fact  is  used  in  the 
above  equation  to  lead  to  the  conclusion  that  kN(co)  <1  V&>  and  that  equation  (2.12)  is 
a necessary  and  sufficient  condition  for  robust  stability. 

It  should  be  noted  that  the  Nyquist  robust  stability  margin  must  be  less  than  1 at 
all  frequencies  and  amplitudes  to  guarantee  stability;  therefore,  kN(co,a)  must  be 
maximized  over  all  frequencies  and  amplitudes  in  order  to  assess  the  robust  stability.  Let 

kN:-  maxkN(co)  (2.13) 

<D>  0 


Then,  from  Theorem  2.2  it  follows  that  a necessary  and  sufficient  condition  for  the  robust 
stability  of  the  uncertain  nonlinear  closed  loop  is  that  kN  < 1 . 
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Theorem  2.2  states  that  the  scalar  kN(a>)  is  used  to  characterize  the  robust 
stability  of  the  closed-loop  system.  The  computation  of  kN(co)  requires  the  knowledge 
of  the  critical  perturbation  radius  defined  in  (2.8).  The  challenging  task  in  this  problem  is 
to  calculate  the  critical  perturbation  radius. 

When  vc(a> ) is  convex,  equation  (2.7)  indicates  that  pc(<x>)  is  the  distance  from 
the  nominal  point  g0(ja>)  and  the  unique  point  where  the  critical  line  intersects  the 
boundary  of  the  uncertain  value  set  u(co) . However,  when  vc((o)  is  non-convex  there 
are  multiple  points  where  the  critical  line  intersects  v(a>) . Then  equation  (2.8)  states  that 
pc(co)  is  a function  of  the  distance  between  g0(jco)  and  the  boundary-intersection  point 
that  is  closest  to  the  critical  loci.  Because  the  convexity  of  uc(a>)  may  not  be  known  a 
priori,  the  general  definition  of  the  critical  perturbation  radius  (2.8)  allows  the  critical 
direction  theory  to  be  applied  to  a more  general  class  of  uncertain  systems,  especially  the 
case  of  real  affine  uncertain  systems. 

2.4.  Special  Cases  Under  Consideration 
2.4.1  Derivations  Using  Circular  Value  Sets 

For  the  circular  value  sets  considered,  the  critical  radius  pc(a>,a)  required  in 

(2.10)  or  (2.11)  for  the  calculation  of  kN{co,a)  is  explicitly  given  by  (2.2),  namely 
pcico,a)  = r(<y)  is  a constant  at  each  frequency  and  is  independent  of  the  amplitude  a . 
Consequently,  the  optimization  (2.10)  or  (2.11)  required  to  determine  the  Nyquist  robust 
stability  margin  can  be  carried  out  in  an  efficient  fashion  by  minimizing  the  distance 


d(a>,a  ) = 


n{a  ) 


+ goUa) 


(2.14) 


from  the  nominal  point  to  the  critical  loci.  Define  now  the  minimizing  amplitude 
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a (w)  = argmin 

a>  0 


n(a ) 


+ g0U®) 


A value  for  kN(co)  can  then  be  obtained  from  (2.10)  or  (2.11)  at  each  frequency  if 
a (jo)  is  known.  Theorem  2.3  illustrated  below  provides  details  of  the  procedure. 

Let  a circular  uncertain  value  set  described  by  a radius  r{co)  and  let  d(oj,a  ) 
denote  the  minimum  distance.  Then  the  Nyquist  robust  stability  margin  is  given  by 


k^i®) ~ 


r(,<o) 

♦ 

d(co,a  ) 


(2.15) 


where  the  minimizing  amplitude  is  given  by 

a*  = min  {d(co  ,a\),  d(co  ,a2),  - ■ ■ , d(a>  ,ak)} 

and  where  the  set  j a*  ,a*2,---,a*k}  represents  all  roots  of  the  equation 


a + f 


(nr(a)2  +«,-(fl)2) 


= 0 


(2.16) 


where 

a(jo,a ) = (~nr (a)2  gor (Jet))  + 2 nr (a)nt ( a)goi (jo)  + nt (a)2  gor  (jeo)  - nr (a))  ^ - (2-17) 

/3(ct), a)  = («,. (a)2  goi (. jeo ) - 2 nr (a)nt ( a)gor (jeo)  - nr(a )2 goi (jeo)  - n, (a))  (2.18) 

Proof.  The  derivative  of  (2.14)  must  be  calculated  to  find  the  minimizing  amplitude 
a* (a>) . First,  rewrite  equation  (2.14)  as 


d(a>,a)  = 


1 


: + 


— - I 1/2 

8°U®)  + goU°)\  + go(jo))go(jco)  i (2.19) 


n(a)n(a ) n(a)  n(a) 
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The  derivatives  of  each  term  inside  the  braces  in  (2.19)  are  determined  as  follows.  Split 
the  nominal  point  and  describing  function  into  their  respective  real  and  imaginary  parts, 
i.  e. , 

go  0®)  = gor  (ja>)  + Jgoi  O'®) 
go0'®)  = SorO®)  - jgoi  O'®) 


and 


fi(a)  = fir(a)  + ./fii(a) 


n(a)  = nr(a)-jni(a) 


Now,  the  derivative  of 


1 


n(a)n(a ) 


can  be  written  as 


A 

da 


r \ -2 

1 


f c/«,.(a)  _ / A dni(a)'' 


|«(a)| 


(«r(a)2  +n1(a)2) 


(2.20) 


The  terms  So^a>\+  S°^JC01  in  equation  (2.19)  are  rearranged  and  simplified  to  yield 
n(a)  n(a) 


gn(jco)  «n(i(o)  2(/ir(fl)gorOffl)-n/(fl)goiO'®)) 

— = 1 — — i . . n 


(2.21) 


n(a)  n(a ) nr(a)  + «,(a)" 

Using  the  quotient  rule,  the  derivative  of  (2.21)  with  respect  to  the  amplitude  can  be 


written  as 


d 

da 


g0U ®)  ! g„0‘®) 
n(a)  «(a) 


SorO'®) 


dnr(a ) 
da 


goi(ja) 


dn,(a) 

da 


(nr(a)2 +ni(af) 


, / , sdnJa)  , .</«.■  (a)  ^ 

4(nr(fl)gorO'®)-”/(fl)S'oiO®))  — 7 + w/(a) 


c/a 


c/a 


2 

(«r(a)2+n,(a)2) 


(2.22) 
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Finally,  the  term  g0(jo))g0{jco)  is  not  dependent  upon  the  amplitude  a , making  its 
derivative  identically  zero.  The  derivative  of  (2.19)  is  readily  obtained  using  equations 
(2.20)  and  (2.22)  and  carrying  out  straightforward  algebraic  manipulations,  resulting  in 

(2.23) 


d ( u \\  1 a + P 

— (J(«,«))  = — -iTj-  772 

da  d(co,a)  ^nr(a)2  +tij(a)  j 


where  a and  /?  are  respectively  given  by  (2.17)  and  (2.18).  Setting  (2.23)  equal  to  zero 
and  recognizing  that  d(co,a)  & 0 by  the  assumption  of  nominal  stability  yields  equation 


(2.16). 

The  robustness  analysis  for  non-circular  templates  is  challenging  because  both  the 
critical  perturbation  radius  and  the  distance  from  the  nominal  point  to  the  critical  loci 
change  as  a function  of  frequency  and  amplitude.  To  simplify  the  calculation,  it  is 
common  to  circumscribe  an  arbitrary  shaped  template  with  an  appropriate  circle 
(Bhattacharya,  1993),  making  the  uncertainty  value  set  u(co)  become  a circle  with  radius 
r{co) . Therefore  at  each  frequency,  the  critical  perturbation  radius  becomes  a constant 
equal  t o t he  r adius  o f t he  c ircle  and  i ndependent  o f t he  a mplitude.  T he  a nalysis  t hen 
proceeds  as  in  the  circular  case.  The  resulting  condition  kN  < 1 is  now  sufficient  only, 
given  that  the  original  value  set  has  been  substituted  for  a . 

* 

Theorem  2.3  can  be  used  to  derive  analytical  expressions  for  the  amplitude  a (co) 
that  solve  the  equation  d(d(a),a))/  da  = 0 for  several  nonlinear  operators  of  interest. 

Table  2.1  shows  four  nonlinear  maps  reported  by  Mohler  (1991),  which  afford  analytical 
treatment.  The  first  column  of  Table  2.1  shows  in  graphical  form  the  nonlinear  operator 
considered.  The  parameter  m represents  the  point  where  the  map  / (e)  intercepts  the 
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ordinate,  a is  the  magnitude  of  the  map  when  it  intercepts  the  positive  abscissa,  and  k 
represents  the  slope  of  the  map.  The  next  two  columns  show  the  real  and  imaginary  parts 
for  describing  function  of  the  map,  and  finally  the  last  column  gives  an  analytical 
expression  for  the  candidate  amplitudes  that  must  be  considered  as  possible  minimizers 
of  (2.14).  For  example,  the  third  row  shows  that  there  are  two  amplitudes  that  set  the 

* -Am  g0„(a>)  , *,  \ 

derivative  (2.16)  equal  to  zero,  namely  ax(co)  = — and  a2{co)-*<n. 

n 1 + kgor(eo) 

Hence,  the  minimizer  is  to  be  selected  from  these  two  candidate  amplitudes. 

The  nonlinear  operator  for  saturation  is  shown  below.  This  nonlinear  element  can 
be  approximated  by  the  following  describing  function 

n(a)=k  for^Sl  (2.24) 

n(a)  = kNu  else 


where 


Figure  2.5.  Graph  for  saturation  illustrating  the  slope  k and  endpoints  d. 
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Table  2.1.  Describing  functions  and  analytical  expressions  for  the 
amplitudes  that  solve  the  equation  d(d(co,a))  I da  = 0 . 


f(e) 

nr(a) 

ni(a) 

Candidate  minimizing  amplitudes 

m . 

4 m 

0 

\ -4 mgor{co) 

ax  ( co ) - 

71 

-m 

na 

m 

, 4 m 

k H 

0 

* N - 4m  gnr{co ) *.  . 

. l , a2(co)- co 

n l + *gor(®) 

y 

-m 

na 

m 

T 

1 

k 

4m 

na 

•,  \ 8 mg0i{(0)  * %mgoi(co) 

a\(G))-  „ , x a r \ ’ 

n6x(a>)  n02(co) 

:|e 

(<w)  - co 

6,1(«)  = 2^gor(«)  + l- 

t 

j 

^k2  (glr  (®)  + gli  (®))  + 4%or  (®)  + 1 

02(o)  = 2&gor(»  + l + 

^4£2  (g2,.  (ffl)  + go«  (®))  + 4^or  (®)  + 1 

The  approach  used  earlier  cannot  yield  a simple  analytic  expression  because  ot 
the  complex  nature  of  the  describing  function.  However,  an  analytical  expression  can  be 
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obtained  analyzing  the  describing  function.  The  analytic  expressions  for  the  minimum 
distance  from  the  nominal  point  to  the  critical  loci  are  given  as  follows 


d*x  («)  = 


f 

gli(®)  + g2or(°})  + 

V 


2ng0r(co)  | 
2.3768& 


( 


n 


V 


2.3768A: 


2\ 


1/2 


<£(®)=koi(®)| 
^3(®)  = ko(®)| 


(2.25) 


Hence,  the  minimizing  distance  <7  (ft))  is  to  be  taken  from  one  of  the  three  candidate 
distances. 

Proof.  The  describing  function  given  by  (2.24)  maps  to  a straight  line  segment  on 
the  Nyquist  plane  as  shown  in  Figure  2.6. 


Im 


Re 


n 


2.3768A: 

Figure  2.6.  Nyquist  map  of  the  saturation  describing  function.  Note  that  the 
endpoint  depends  only  on  the  slope  of  the  saturation  line. 

There  are  now  three  regions  where  the  nominal  point  g0(a>)  may  lie:  (1)  to  the  left  of  the 
line  segment,  (2)  above  or  below  the  line  segment,  or  (3)  to  the  right  of  the  line  segment 
and  we  will  treat  each  case  in  turn.  If  g0(ft>)  lies  to  the  left  of  the  line  segment,  then  the 


shortest  distance  will  be  from  g0(<w)  to  the  endpoint  - 


n 


2.3768/c 


which  is  given  by 


d\  (©)  = | g0i  («)|2  + 8 or  (®)  + 


n 


\ 


2.3768A: 


which  can  be  simplified  to  yield  d\  (oj)  . When  the  nominal  point  lies  directly  above  oi 
below  the  line  segment,  region  (2),  the  shortest  distance  will  simply  be  the  magnitude  of 
the  imaginary  part  of  the  nominal  point  goi(a> ) • The  last  area  to  analyze  is  when  g0(a>) 
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lies  to  the  right  of  the  line  segment.  The  shortest  distance  will  then  be  the  distance  from 
the  origin  to  the  nominal  point  g0(a> ) . Thus  the  minimizing  distance  is  given  by  the 

* 

magnitude  of  the  nominal  point  illustrated  by  d 3 ( co ) . 

2.4.2  Real  Affine  Parametric  Uncertainty 

The  generalized  critical  direction  theory  developed  in  Section  2.3  is  specialized  to 

systems  with  real  affine  parametric  uncertainties  of  the  form 

p 

n0(s)  + J]qini(s) 

g(s,  q)  = y , q e Q (2-26) 

do(s)  + Yjqidi(s) 

i= 1 

where 

l 

«o0):=  Zwo  ksk 

k= 0 
and 

m 

do(s)\=  Yjdoksk 
k= 0 

are  known  nominal  polynomials, 

t 

>h(s) = Yjn^k 

k= 0 

and 

m 

d,(s)  = 

k=0 

are  known  perturbation  polynomials,  and  q = [<?i  ^2  •••  e^P  i s a v ector  o f real 

perturbation  parameters  belonging  to  the  bounded  rectangular  polytope 
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Q = {qeRp 


Qi  —Qi—Qi  ’ 2 1,2,. 


(2.27) 


where  <7,:  and  q* , i = 1,  2,  ...,  p are  finite  real  bounds.  Equations  (2.26)  and  (2.27) 
define  a class  of  finite  dimensional,  linear,  time-invariant,  real  plants  with  affine 


uncertainties. 

The  first  step  in  determining  the  generalized  critical  perturbation  radius  for  the 
uncertain  system  (2.26)  - (2.27)  is  to  determine  whether  the  critical  loci  -\/n(a)  belongs 


to  the  uncertainty  set  o(co)  . 

The  affine  uncertain  system  (2.26)  can  be  written  in  vector-matrix  form 


Zr 


1 5 


/-1  5' 


g(s,q)  = 


”00 

«01 

.”0/ 


+ 


n\0  n20 
«11  n2\ 

«1  / n2l 


np0 

np\ 

npl 


q\ 

<12 


1 5 


sm~l  sm 


^00 

^01 


d 


0 m 


+ 


<^10  ^20 
dn  d2\ 


d 


1 m u2m 


lp  0 
1p\ 


1 pm 


<i\ 

d2 


sI(no  +N^) 


SI( 


(2.28) 


;(d0+Dpq) 

where  s n and  are  vectors  of  lengths  i + 1 and  m + 1 , containing  powers  of  the  Laplace 
variable  s,  and  where  n0  eR£+1,  do  e^"+1>  ^ ^ P anc*  ^ are 


constant  vectors  and  matrices  that  represent  the  structure  of  the  affine  parametric 
uncertainty.  The  value  set  (2.28)  is  evaluated  at  a frequency  co  by  substituting  5 = jco  to 


yield 
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g(ja>,  q) 


smI 

[no  ,/t+Npvq) 

\ + j*n,l\ 

no,/+NP(/q) 

sd,R\ 

(d0,/f  +Dp,/f(l)  + 7sJ,/ 

(do,/ +D/>,/q) 

(2.29) 


where 

*M=[1 

r}  ("  7/21+1  T 3 5 7 

•••  eR1  1 ,snI=  co  -co  co  -co 

./  eRr(M/2i 

= 1 -«2  CO4  -CO6 

T.rm/21+l  T 3 5 7 

•••  gR1  1 ,s  dI-  CO  -CO  CO  -CO 

eRfK0/2i 

and 


”00 

n0,/f  - 

1 

a 

...  0 
to 

sR^1, 

Npv  = 

”01 

n0,7  = 

”03 

£Rr('+i>'ji, 

Np,/  = 

O 

O 

^3 

- 

d0,7?  = 

^02 

Dp,R  = 

^01 

eRl(m+»/2l 

d0,7  = 

^03 

D,'  = 

”10  ”20 
”12  ”22 


”11  ”21 
”13  ”23 


^10  ^20 
^12  ^22 


J21 

J,3  d22 


np0 

np2 


np  1 

”p3 


'pO 

h>2 


gR 


(R/2>1)> 


eR 


f(^l)/2>p 


eR 


(|"m/2~|+l)x/? 


pi 


P3 


eR 


|"(m+l)/2"|xp 


where  |~ •]  represents  the  greatest  integer  function. 

Now  consider  a point  on  the  critical  loci  w-wR+  jwj  g-1  / n(a) . It  can  be 
shown  that  w g u(co)  if  and  only  if  there  exists  a vector  q eQ  such  that  g{jco,q)  - w. 
Using  (2.29)  this  requirement  is  equal  to  obtaining  a vector  q eQ  that  solves  the 


equation 
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s Tn,R  (n  0,R  + Np,*q)  + jsTnJ  (n  0>/  + N pJq) 

— - ! — r 7jT~) T = WR+JWj 

sd,ll(A0,R  +Dp,«q)  + >rf,/(dO,/  +Dp,/tl) 


(2.30) 


This  can  be  characterized  by  a linear  equality  problem  shown  by  Theorem  2.4: 

Let  w 6 — 1 / n(a)  point  with  finite  magnitude  on  the  critical  loci.  Then 

w e o(a>)  if  and  only  if  there  exists  a feasible  solution  q g Rp  to  the  linear  equality 


A(w)  q = b(w) 

subject  to  the  linear-inequality  constraint 


(2.31) 


1 

0 

0 •• 

o' 

4i+ 

-1 

0 

0 •• 

0 

-q\ 

0 

1 

0 •• 

0 

42 

0 

-1 

0 •• 

0 

q^ 

-42 

0 

0 

0 •• 

1 

4; 

0 

0 

0 •• 

• -1 

rq~p. 

(2.32) 


where 


A(w):  = 


sn,R^p,R  -wRsd,R^p,R+wIsd,I^pJ 

S IjKpJ  - wRSd,IDp,l  ~ WISd,R^p,R 


b(w):  = 


'sl,Rn0,R  +wRsd,R^0,R~wISd  ,1^0,1 


gR 


G R 


2x  p 


(2.33) 


(2.34) 


~sn,In0J  +wRsd,IA0,I  + wISd,Rd0,R 

Proof  A finite  point  w e-l/n(a)  belongs  to  v(co)  if  and  only  if  there  exists  a vector 
q e Q that  satisfies  equation  (2.30).  Since  the  denominator  of  the  left  hand  side  of  (2.30) 
is  non-zero  because  of  the  finite  magnitude  of  w , the  equality  can  be  rearranged  and  the 
real  and  imaginary  parts  separated  to  yield 
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sl,R^p,R  ~ wRsd,R^p,R  + wIsd,I®p,I 
sn,I^p,I  ~ wRsd,I D p,I  ~ wIsd,R®p,R 


~sn,RnO,R  + wRsd,RdO,R  " wIsd,I^O,I 
_s«,/n0,/  +wRsd,IdOJ  + wIsd,RA0,R 


which  is  equal  to  (2.31)  after  the  definitions  given  in  (2.33)  and  (2.34)  are  substituted. 
From  the  definition  of  Q'm  (2.27),  the  constraint  that  qeQ  can  be  described  by  the 
linear  inequality  (2.32).  Therefore,  w e-l  / n(a)  is  an  element  of  u(a>)  if  and  only  if 
there  exists  a feasible  solution  to  the  linear  program  defined  by  equations  (2.31)  and 
(2.32). 


The  computation  of  the  Nyquist  robust  stability  margin  kN(co)  for  affine  systems 
(2.26)  requires  that  the  critical  perturbation  radius  pc(co)  be  calculated  first.  As 
illustrated  by  (2.8)  and  (2.9),  this  requires  determining  the  set  j3c(co).  Once  all  the 
elements  of  /3c(co ) have  been  determined,  it  is  straightforward  to  calculate  pc{co)  from 
formulas.  For  the  case  of  affine-uncertain  systems  of  the  form  given  by  (2.26),  Baab  et 
al.  (2001)  have  identified  a two-step  strategy  that  effectively  determines  the  critical 
boundary-intersections  set  Pc(co) . For  brevity  it  has  been  omitted  here. 

2.5.  Examples 

2.5.1  Circular  Value  Set  with  Simple  Describing  Function 

Consider  the  open-loop  unstable  system 


1 

-5  + 10 


with  an  associated  unstructured  (circular)  uncertainty  with  radius 


1 

-j(D  + 10 


r{a)  = 
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and  a nonlinear  operator  f ( e ) of  the  form  given  in  the  last  row  of  Table  2.1  with  k = 2 
and  m = 1 . The  describing  function  in  this  case  is 

n(a)  = 2 +—j 
na 

It  is  straightforward  to  show  that  the  nominal  closed-loop  system  (i. e. , when 
r(co)  - 0 \fa> ) is  stable. 

Figure  2.7  shows  the  nominal  point  ga(jco)  at  a frequency  co  - 10 , along  with  a 
circular  disk  representing  the  uncertain  value  set  at  this  frequency.  At  this  frequency, 
pc(a>,a)  = r(co)  = 0.0707  . The  figure  also  shows  the  critical  loci  for  amplitudes  ranging 
from  0.0001  to  0.5. 


i l i \ I 1— 
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Figure  2.7 . Nominal  point,  uncertainty  value  set,  and  critical  loci  for  the 
example. 

From  Table  2.1  (last  row,  last  column),  the  three  candidate  amplitudes  that  must  be 


considered  are 
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a\ (a)  = 


8 mgoi(eo) 
n6x(co) 


-7.692 


a*2(co)  = 


8 mgoi(aj) 
k62{o}) 


0.0527 


and 


* 

a3(a>)  = oo 

Noting  that  d(a>,  - 7.692)  does  not  exist  because  the  amplitude  must  be  strictly  positive, 
d(a>,  0.0527)  = 0.0541 , and  d(co,  oo)  = 0.5523 , it  follows  that  the  minimum  distance 

sought  is  achieved  at  a*  = 0.0527.  From  equation  (2.15)  it  follows  that 
kN(co)  = 1.306  > 1 . It  can  be  concluded  from  Theorem  2.2  that  the  uncertain  nonlinear 
closed  loop  system  is  not  robustly  stable  with  respect  to  the  given  uncertainty. 

Further  insight  is  gained  from  Figure  2.8,  which  shows  a plot  of  the  argument  of 
the  right-hand  side  of  (2.10)  as  a function  of  the  amplitude  a.  The  infimum  value 
* 0.1280  occurs  as  a — > oo  and  the  maximum  value  = 1.306  occurs  at  a = 0.0527 . The 
numerical  results  are  consistent  with  the  numerical  root-analysis  discussed  above. 
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Figure  2.8.  The  Nyquist  robust  stability  margin  plotted  versus  the  amplitude 
for  the  example. 

2.5.2  Circular  Value  Set  with  Saturation  Describing  Function 

Consider  the  open-loop  stable  system 


, \ 25 

= -2 7 

5-5  + 1 


with  an  associated  unstructured  (circular)  uncertainty  with  radius 


(. ja>) 2 -jto  + l 

and  a nonlinear  operator  f(e)  of  the  form  given  in  Figure  2.5  with  k = 3 and  d = 2. 
The  describing  function  in  this  case  is  given  by  equation  (2.24).  It  is  straightforward  to 
show  that  the  nominal  closed-loop  system  is  unstable. 

Figure  2.9  shows  the  nominal  point  g0(j<o ) at  a frequency  a>  = 0.8 , along  with  a 
circular  disk  representing  the  uncertain  value  set  at  this  frequency.  At  this  frequency, 


r(co)  = 
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pc(co,a)  = r(a> ) = 0.6079 . The  figure  also  shows  the  critical  loci  for  amplitudes  ranging 
from  0.01  to  10,000. 


Figure  2.9.  Nominal  point,  uncertainty  value  set,  and  critical  loci  for  the 
saturation  describing  function. 

Because  the  nominal  point  is  to  the  left  of  the  critical  loci,  the  smallest  distance  given  by 
(2.25)  is 


d\  (®)  = 


Soi(C0)  + Sor(C0)  + 


2xgor(v) 

2.3768A: 


+ 


n 


2.3768£ 


2^\ 


1/2 


= 1.4335 


From  equation  (2.15)  it  follows  that  kN  (0.8)  = 0.4241  < 1 so  the  system  is  robustly 
stabilizing  at  the  given  frequency. 

However,  Theorem  2.2  states  that  kN(co)  <1  \f  co  in  order  for  the  system  to  be 


robustly  stable  for  all  frequencies.  Figure  2.10  shows  the  values  of  kN(co ) calculated  for 


a sequence  of  points  spaced  in  the  range  [0.001,  10].  The  maximum  value  of  the  Nyquist 
robust  stability  margin  is  kN(a>)  - 0.4241  which  occurs  at  a frequency  <y  = 0.8.  From 
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Figure  2.10  it  is  readily  concluded  that  kN(a>)  < 1 therefore  the  closed  loop  is  robustly 
stable  for  the  given  uncertainty. 


Figure  2.10.  The  Nyquist  robust  stability  margin  plotted  versus  the 
frequency  for  the  saturation  example. 

2.5.3  Affine  Uncertainty  Structure  with  Describing  Function 

Consider  the  system  with  affine  uncertainty  structure 

, (0.6s3  + 2.9s2  + 7.5  s + 15)  + (0.1 2s2  + 0.7s  + 1 )o.  + (0.06s2  + 0.2  s)q2  + (-0.3s  - 1 )q, 

g(s,  cj)  — 

(2s4  + 1 Is3  + 30s2  + 20s  + 0.2)  + (0.5s3  + 2s2  - s)ql  + (-0.5s3  + s2  )q2  + (0.5s3  + s)q2 


where  the  parameter  uncertainty  vector  belongs  to  the  rectangular  polytope 


Q={qeR3 


-10<?,  <10, 


-0.3<q2<03, 


- 0.3  < q2  < 0.3} 


With  this  uncertainty  the  system  now  has  a non-convex  critical  uncertainty  value  set 
vc(a> ) at  the  frequency  (0  - 1.6,  as  can  be  determined  by  inspection  of  the  frame 
g(jco,E(Q))  of  the  value  set  shown  in  Figure  2.11.  The  figure  shows  that  as  one  travels 
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from  the  nominal  point  g0(co)  in  all  directions  spanned  by  the  critical  cone  9c{co)  that 
there  are  s everal  i ntersections  w ith  t he  v alue  s et  b oundary.  A s a c onsequence  o f t his 
geometry,  the  critical  uncertainty  value  set  uc(o))  is  a collection  of  lines  and  the  unions 
of  two  disjointed  line  segments. 


0.1  r 
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Re 

Figure  2.11.  Frame  for  the  uncertainty  set  for  the  non-convex  system  at  the 
frequency  co  = 1.6 . 

Theorem  2.4  is  used  to  determine  the  linear  feasibility  for  every  point  belonging 
to  the  critical  loci-l/n(a) . For  a frequency  a>  - 1.6,  it  can  be  readily  verified  that  the 
problem  is  infeasible  for  all  points  along  the  critical  loci.  From  Theorem  2.4  it  follows 
that  -1  /n(a)<£u(a>)  and  therefore  it  can  be  claimed  that  at  co-  1.6,  the  value  set 
excludes  the  critical  loci. 

The  critical  perturbation  radius  pc(co)  and  Nyquist  robust  stability  margin 
kN(a>)  can  now  be  calculated  for  any  frequency  using  the  method  described  by  Section 
2.4.2  and  Baab  et  al.  (2001).  There  are  an  infinite  amount  of  directions  to  search  inside 
the  critical  cone  9c(co)  and  therefore  only  a finite  amount  of  critical  directions  can  be 
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searched.  The  Nyquist  robust  stability  margin  is  calculated  for  each  critical  direction  and 
the  largest  kN(co)  that  results  is  employed.  Invoking  equation  (2.8)  it  is  readily 
determined  that  at  the  frequency  co  = 1.6 

M1.6)  = 0-7698  <1 

which  is  a result  that  is  consistent  with  the  conclusion  reached  using  Theorem  2.4  that  the 
value  set  does  not  include  the  critical  loci  at  this  frequency. 

Figure  2.12  illustrates  the  values  of  kN(a> ) calculated  for  a sequence  of  200 
frequency  points  equally  spaced  in  a logarithmic  scale  in  the  range  [0.001,  10]. 

1.2  - | 
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( O 

Figure  2.12.  Nyquist  robust  stability  margin  kN(co)  for  the  real  affine 
parametric  uncertainty  system 

From  Figure  2.12  it  is  readily  concluded  that  kN(co)  > 1 making  the  closed  loop  unstable 
with  respect  to  the  given  uncertainty. 

2.6.  Conclusions  and  Future  Work 

This  paper  extends  previous  work  on  robustness  theory  for  linear  systems  to  the 
case  of  nonlinearities  that  can  be  captured  through  the  describing  function  method.  A 
new  definition  of  the  Nyquist  robust  stability  margin  for  both  convex  and  non-convex 
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systems  takes  into  account  a critical  loci  via  the  describing  function  approach  is 
introduced.  Circular  value  sets  are  studied  and  analytic  expressions  for  candidate 
amplitudes  for  minimizing  the  Nyquist  robust  stability  margin  are  derived  for  several 
nonlinear  examples  including  the  case  of  valve  saturation.  Real  affine  parametric 
uncertain  systems  were  analyzed  via  a feasibility  problem  and  arc  intersections  similar  to 
Baab  et  al.  (2001).  Finally,  three  examples  were  given  that  illustrated  the  analytic 
expressions  obtained  and  the  non-convex  definition  of  the  critical  perturbation  radius. 


CHAPTER  3 

SYSTEMATIC  ALGORITHMS  FOR  CHARACTERIZING  THE  ROBUSTNES  OF 

PREDICTIVE  CONTROLLERS 

3.1.  Introduction 

Linear  models  of  plants  used  for  designing  predictive  controllers  are  affected  by 
different  degrees  of  uncertainty  that  normally  appear  as  errors  in  the  coefficients  to  the 
numerator  and  denominator  polynomials.  Such  real  parametric  uncertainties,  although  of 
unknown  value,  are  often  known  lie  within  a domain.  Representative  uncertainty-domain 
types  include  the  cases  where  the  domain  can  be  described  as  a hyperellipsoid  or  as  a 
hyperbox  (i.e.,  a rectangular  polytope). 

More  s pecifically,  E llipsoidal  u ncertainties  a re  d efined  a s a v ector  o f u ncertain 
plant  parameters  that  are  constrained  to  lie  inside  a known  ellipsoid.  These  descriptions 
are  common  when  using  regression  techniques  to  fit  experimental  data  to  model 
parameters  (Belfonte  et  al.,  1990),  and  are  often  encountered  in  the  context  of  predictive 
control  design.  These  regression  techniques  produce  a nominal  value  for  the  parameters 
and  a variance-covariance  matrix  that  is  used  to  specify  an  ellipsoidal  r egion  centered 
around  the  nominal  plant  values.  Ellipsoidal  descriptions  have  the  advantage  that  they 
can  be  quantified  using  systematic  statistical  methods,  as  opposed  to  previous 
unstructured  uncertainty  descriptions  which  are  difficult  to  quantify  (Kouvaritakis  et  al., 
1992;  Hrissagis  et  al.,  1996)  and  often  depict  excessively  conservative  domains  (Zafirou, 
1990;  Genceli  and  Nikolaou,  1993).  Hyperbox  uncertainty  descriptions  appear  when  the 
parameter  identification  procedure  yields  upper  and  lower  bounds  for  each  of  the  nominal 
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parameter  coefficients.  Such  descriptions  can  also  be  obtained  from  the  statistics  of  the 
parameter  identification  procedure,  from  physical  arguments  (for  example,  knowledge 
that  t he  m aximum  e rror  p ossible  i n the  p arameter  i s w ithin  p lus  o r m inus  1 0%  o f t he 
nominal  value). 

The  challenge  for  robust  predictive  control  design  is  to  find  a set  of  tuning 
parameters,  namely  control  weights  and  prediction  horizons,  that  produce  a stable  closed 
loop  for  all  the  plants  belonging  to  a given  real  parametric  uncertainty  description.  This 
chapter  presents  a methodology  for  (i)  robust  stability  analysis,  i.e.,  determining  whether 
a given  predictive  control  design  is  robustly  stabilizing  with  respect  to  the  given  real 
parametric  uncertainty  domain,  and  (ii)  robust  predictive  control  synthesis,  i.e.,  finding  a 
set  of  tuning  parameters  that  ensure  that  the  design  is  robustly  stabilizing  with  respect  to 
the  given  uncertainty  domain.  The  method  proposed  is  based  on  a numerical  approach 
that  allows  addressing  in  a systematic  the  robustness  analysis  and  robust  synthesis 
problems.  A numerical  approach  is  adopted  because,  unlike  previous  robust  design 
methodology  for  predictive  controllers  using  the  Youla  parameter  (Mahon  et  al,  1999), 
the  underlying  optimization  problems  are  not  convex.  The  computational  burden  is 
nevertheless  r easonable  for  off-line  calculations,  and  the  underlying  algorithms  d o n ot 
require  significant  effort  in  developing  software  solutions  because  they  involve  standard 
and  readily  available  tools,  such  as  linear  programming,  for  example.  The  approach 
allows  t he  c ontrol  d esigner  t o e ffectively  i dentify  predictive  c ontrol  t uning  p arameters 
that  improve  the  robust  stability  of  the  loop. 
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3.2.  Formulation  of  the  Problem 
3.2.1  Uncertainty  Description 

Consider  the  nominal  plant  transfer  function 


y(z)  = 


B(z) 

A(z) 


u(z) 


(3-1) 


wth  input  u(z) , output  y(z) , and  with  known  coprime  polynomials 

A(z)  = zn  + an_\zn  * + •••  + «() 
and 

B{z)=bmzm  +bm_xzm~l  +---+b0 

of  order  n and  m,  respectively,  where  n>m.  The  uncertain  process  model  is  given  by 
the  transfer  function 


y(z)  = 


5CZKA5CZ) 
A(z)  + AA(z ) 


(3.2) 


where  the  polynomial  uncertainties 

AA(z)  = San_\zn~^  +San_2Zn  2+---  + SaQ 
AB(z)  = Sbmzm  + 5bm_\zm 

feature  uncertain  real  coefficients  defined  by  the  vector 

<Sp  = [ San_i  San_2  ■■■  Sqq  Sbm  Sbm_\  ■■■  5b 0]7  eiR,i+"'+1  (3.3) 


>n+m  + 1 


that  belongs  to  a known  closed  parametric  domain  A , i.e.,  <5p  <=  A e iK 

The  generality  of  the  approach  is  enhanced  if  the  description  of  the  uncertainty 
domain  can  be  cast  in  the  normalized  form 


■{ 


A = <te<R 


n+m+ 1 


**} 


(3.4) 
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where  ||•||  represents  the  vector  p-norm.  The  choice  of  an  appropriate  norm  depends  on 

the  uncertainty  description.  One  case  of  interest  is  where  the  uncertainty  domain  is  the 
ellipsoid 


which 
Q eiR 


Ae  = {<Sp  < l} 


is  characterized  by  a known  symmetric  positive-definite 
(n+m+l)x(n+m+l)  Producing  the  coordinate  transformation 


(3.5) 

matrix 


(3.6) 


allows  rewriting  (3.5)  as 


Ae 


G^+m+1 


(3.7) 


which  is  in  the  normalized  form  (3.4). 

A second  case  considered  is  where  the  uncertainty  description  is  given  as  the 
symmetric  rectangular  polytope 

A^  = g93”+,”+1  I — b < < b}  (3.8) 

which  is  fully  characterized  by  the  vector  b eMn+m+l  that  contains  positive  elements. 
The  inequalities  in  (3.8)  are  interpreted  to  apply  in  an  element-wise  fashion.  Consider  the 
coordinate  transformation 

th  = diag(b)~^Sp  (3.9) 

where  dicig( b)  e <j^(w+/m+1)x(«+/m+1)  ^ diagonal  matrix  with  diagonal  elements 

b\  ,b2 , ■ ■ ■ , bn + m + 1 . Then  the  hyperbox  domain  (3 . 8)  can  be  rewritten  as 

At={t4€«"+”’+1|||tt|si}  (3.10) 
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which  is  of  the  normalized  form  (3.4).  The  transformations  (3.6)  and  (3.9)  can  be 
represented  in  the  generic  form 

t = D%  (3-11) 

where  D = Q1/2  and  t = te  in  the  case  where  A = Ae,  and  t = tb  and  D = diag{ b)  in 
the  case  where  A = . 

3.2.2  Robustness  of  Predictive  Control  Designs 

A controller  is  said  to  be  a predictive  controller  if  at  every  sampling  instant  t it 

minimizes  the  performance  functional 

N,  N . 

J(t)  = Y,[r(t  + i)  - y(t  + i\t)]2  +/lX[Aw(t  + 01  (3-12) 

i=l  i=0 

where  {/*(/  + z')}  is  a sequence  of  future  setpoints,  {y(t  + /|0}  is  the  series  of  future 
values  of  the  output  predicted  at  time  t,  and  {Aw(/+z)}  is  the  sequence  of  future 
controller  inputs.  Tuning  is  achieved  by  selecting  appropriate  values  for  the  control 
weighting  parameter  X , the  prediction  horizon  Ny , and  the  control  horizon  N u . 

Crisalle  et  al.  (1989)  have  shown  that  the  predictive  controller  that  minimizes 
(3.12)  can  be  written  in  a transfer-  function  form  (see  also  Mahon  et  al.,  1999;  Hrissagis 
et  al.,  1996).  More  specifically,  in  the  absence  of  uncertainty  the  plant  (3.2)  reduces  to 
the  nominal  plant  (3.1),  and  the  control  law  that  minimizes  (3.12)  is  of  the  form 

^^u(z)  = T(z)r(z)-^^-y(z)  (3.13) 

zn  zn 

where  R{z) , S(z) , and  T{z)  are  controller  polynomials  given  by 

n- 1 


R(z)  = zn  + rn_\z 


+ — l-  tq 


(3.14) 
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S{z)  = snzn  +5„_iz"  1+---  + 50 

(3.15) 

T{z)  = tN  ZN’  +tN  -\ZN’  + • • • + tyz 

y y 

(3.16) 

The  controller  polynomials  R(z) , S(z) , and  T(z)  are  functions  of  the  tuning  parameters 
A,  Ny,  and  Nu,  and  of  the  nominal  polynomials  A(z)  and  B(z) . Furthermore,  these 

polynomials  fulfill  the  requirements  that  R(  1)  = 0 and  7\1)  = S(l),  which  guarantees  the 
presence  of  an  integrator  in  the  predictive  control  law  (Hrissagis  et  al.,  1996),  and 
ensures  that  the  closed  loop  system  realizes  zero  offset  in  servo-response  problems  when 
the  1 oop  i s s table.  T he  resulting  p redictive  c ontrol  s cheme  i s i llustrated  i n F igure  3 . 1 
where  the  controller  is  deployed  in  feedback  with  the  uncertain  plant.  Details  for 
determining  the  controller  polynomials  (3.14)  - (3.16)  can  be  found  in  Crisalle  et  al. 
(1989),  Mahon  et  al.  (1999),  and  Hrissagis  et  al.  (1996). 

Carrying  out  straightforward  block-diagram  algebra  in  Figure  3.1  for  the  nominal 
case  where  there  is  no  uncertainty  ( i.e .,  A A(z)  = A B(z)  = 0)  yields  the  characteristic 
polynomial  of  the  nominal  closed  loop  system 

A*Q{z)  = A{z)R{z)  + B{z)S{z)  (3.17) 

which  can  be  used  to  assess  the  stability  of  the  closed  loop. 

Let  S denote  the  set  of  all  Schur  polynomials  {i.e.,  the  set  of  polynomials  whose 
roots  lie  strictly  inside  the  unit  circle).  The  closed  loop  with  predictive  control  shown  in 

Figure  3.1  is  said  to  be  nominally  stable  if,4o(z)  eS . When  uncertainty  is  included,  the 
uncertain  characteristic  polynomial  is 


A*(z)  = (A(z)  + AA{z))R{z)  + [B{z)  + AB(z))S(z) 


(3.18) 
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and  the  closed  loop  is  said  to  be  robustly  stable  if  A*  (z)  € S for  all  £p  e A . A predictive 
controller  is  robustly  stabilizing  if  it  ensures  the  nominal  and  robust  stability  of  the  closed 
loop. 


Figure  3.1.  Structure  of  a predictive  controller  characterized  by  the  operator 
polynomials  R(z) , S(z ) , and  T(z ) . The  loop  variables  include  the  set 

point  r,  the  input  u,  and  the  exogenous  disturbance  d. 

A scalar  measure  of  the  robust  stability  of  the  closed  loop  is  given  by  the 

parametric  stability  margin 


a - i min  t 
It  eR p 


A (z)  es 


(3.19) 


where ||•||  represents  the  vector  /7-norm  and  where  t — te  for  the  case  A Ae,  and 


t = tb  for  the  case  A = Ab . The  parameter  a is  interpreted  as  the  size  of  the  smallest 
destabilizing  parameter  perturbation.  Note  that  if  a > 1 , the  size  of  the  smallest 
destabilizing  parametric  uncertainty  is  ||t||^  = a > 1 . Given  that  the  set  A includes  only 


uncertain  parameter  vectors  whose  size  is  ||t||  <1,  it  follows  that  no  element  of  A is 

destabilizing.  Hence,  when  a > 1 the  predictive  control  loop  is  robustly  stable  with 
respect  to  all  uncertainties  in  the  set  A . Conversely,  when  a < 1 , there  is  at  least  one 
element  of  A that  is  destabilizing,  and  the  predictive  control  loop  is  not  robustly  stable 
with  respect  to  A.  The  choice  of  /7-norm  in  (3.19)  depends  on  the  structure  of  the 


42 


parametric  uncertainty  domain  A . Obviously,  the  2-norm  is  used  to  analyze  the  case  of 
the  ellipsoidal  uncertainty  domain  (3.7)  and  the  oo-norm  to  treat  the  hyperbox  domain 
(3.10). 

3.3.  Computation  of  the  Parametric  Stability  Margin 

The  focus  for  the  robust  stability  analysis  is  the  characteristic  equation 

A*  (z)  = 0 , or  equivalently 

A4(z)*(z)  + Afl(z)S(s)  = -4J(z)  (3-20) 

which  is  obtained  after  setting  (3.18)  equal  to  zero,  rearranging  terms,  and  using  the 
definition  of  Aq(z)  given  in  (3.17).  The  analysis  of  robust  stability  requires  the 
characterization  of  the  parametric  uncertainty  set  that  places  the  closed-loop  poles  on  the 
unit  circle  z = eJC0,  co  e[0,;r],  namely,  the  set  that  satisfies  A*(eJCJ)  = 0 for  some 
values  of  frequency.  Using  a vector-matrix  representation  of  equation  (3.20),  the  key  to 
the  robust  stability  analysis  is  to  identify  vectors  dp  that  satisfy 

cT(z)Sp  = -Aq(z)  (3.21) 


where 


c(z):= 


zn~lR(z)  zn~2R(z)-R(z)  zmS(z)  zm~lS(z)  ■■■  S(z) 


and  where  z = eJ0) . Furthermore,  using  the  normalizing  transformation  (3.1 1),  equation 
(3.21)  can  be  written  as 

c7’(z)Dt  = -^o(z) 

and  the  calculation  of  the  parametric  robust-stability  margin  (3.19)  then  reduces  to 
solving  the  frequency-dependent  program 
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a(v)  = \mm\\t\\p 


cV'^Dt  = -A^03)  c V®  )Dt  = ~Ao(ejw)\  (3  22) 


tesr 


and  then  setting 


a 


( min  I a(co)  1 
l<ae[0,  J 


3.3.1  Ellipsoidal  Uncertainty 

When  the  uncertainty  is  constrained  to  lie  within  the  ellipsoid  (3.7),  equation 
(3.21)  can  be  expressed  as 

cj(z)te  =-Aq(z) 


where  cj (z)  = c^(z)Q1/2  , and  where  te  is  defined  in  (3.6).  Adopting  the  2-norm  in  the 
definition  of  the  parametric  stability  margin  leads  to  an  optimization  problem  of  the  form 
(3.22)  as  follows: 


a(w)=  minlM2 

t e'Ji "" " 


JCQ 


(3.23) 


The  solution  to  the  standard  linear-quadratic  problem  (3.23)  is  given  by  analytical 


expression 


-Ao(e^) 


e,mm  ~ cTe  (eJ®  )ce (e-i® ) 6 


{ejco) 


so  that 


a(co)  = t 


e,min  H2 


Obviously,  the  uncertain  system  is  robustly  stable  if  and  only  if  a - mina(<y)  > 1 . It  is 

co  e[0 ,n\ 


interesting  to  note  that 


a(0)2  =B(\)2</> 


(3.24) 
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where 


<j,= 


1 


(3.25) 


p(m  + 1) 

for  the  case  where  the  domain  (3.5)  is  a sphere  with  radius  p,  and  for  the  case  of  an 
arbitrary  ellipsoid 


2 

m+ 1 

2 

Il^rsj 

+ 

Z <lss,i 

i=i 

2 

i = l 

2 

9 

rm  + 1 m + 1 3 

Z XQ  SS,ik 
\i= 1 k=\ 


(3.26) 


where  the  submatrices  QRRe$\nxn,  e^m+^xn , Qrs  e9?”x^/”+1Z  and 

Qss  e <>R(,n+1)x(m+1)  are  extracted  from  the  partition 


Q:  = 


Q RR 
_Q  SR 


Q RS 
Qss . 


and  the  vectors  q fS  and  q fS  for  i = 1,  2, m+1  are  extracted  from  the  partitions 


n1/2-- 
Qrs  ■- 


RS  RS 

RS  1 

r*l/2. 

r ss  ss 

SS 

qi  Q2 

■<\m+ 1 

. Qss  := 

Qi  Q2  ‘ 

••  q»i+i 

When  the  ellipsoid  is  a sphere  characterized  by  Q = pi , p>  0 , then  Qrs  and 
Qss  - pi,  and  hence  (3.26)  reduces  to  (3.25). 

It  is  interesting  to  note  that  ar(O)  is  constant  and  is  dependent  only  on  the  plant 
zeros.  The  robust  stability  of  the  system  at  the  frequency  co  - 0 , therefore,  may  be 
altered  by  changing  the  zeros  of  the  nominal  plant. 
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3.3.2  Hyperbox  Uncertainty 

Constraining  the  uncertainty  to  lie  within  the  hyperbox  (3.8),  and  adopting  the  oo- 
norm  in  the  definition  of  the  parametric  stability  margin,  it  is  possible  to  proceed  as  in 
Section  3.3.1  to  formulate  the  following  optimization  problem  of  the  form  (3.22): 

a(ffl)  = |minML 

l t,e9T-“ 


csV‘")t6=-/f0V'”)l 


(3.27) 


where 

cl(ei0)):=cr(ej<u)diag(b) 

The  optimization  problem  (3.27)  has  been  formulated  following  the  development 
in  Bhattacharya  (1995),  where  it  is  shown  that  a solution  can  be  computed  effectively 
using  readily  available  linear-programming  algorithms.  The  solution  t£mjn  is  such  that 


cc(cd)  - min 


Obviously,  the  system  is  robustly  stable  if 

a = mina(<y)  > 1 

i o e [0,  n ] 

In  contrast  to  the  ellipsoidal  case,  an  analytical  expression  for  a(0)  cannot  be  readily 
obtained  because,  in  general,  the  problem  does  not  afford  an  analytical  solution. 

3.4.  Robustification  of  Predictive  Controllers 
Two  different  controllers  are  examined:  (i)  a standard  predictive  controller  (SPC) 
design  that  ignores  uncertainty  ( i.e .,  A A(z)  = A B(z)  = 0)  and  that  is  based  on  equations 
(3.13)  - (3.16)  and  (ii)  a robust  predictive  controller  (RPC)  that  stabilizes  the  loop  for  all 
plants  1 ying  w ithin  a given  u ncertainty  s et  A . T he  p erformance  o f e ach  c ontroller  i s 
evaluated  via  simulation  for  both  a servo  and  regulation  test  for  the  closed  loop  uncertain 
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system.  The  robustification  process  consists  of  analyzing  the  effect  of  the  tuning 
parameters  A,  N y , and  Nu  on  the  robust  stability  margin  a . Only  those  tuning 

parameters  that  ensure  that  a > 1 are  accepted  as  robustly-stabilizing  choices,  provided 
that  such  choices  also  ensure  the  nominal  stability  of  the  loop. 

3.4.1  Ellipsoidal  Uncertainty 

Consider  the  stable  nominal  plant  model 

B ^ (3.28) 

A(z)  z2  -0.8  + 0.16 

with  the  ellipsoidal  parametric  uncertainty  description  (3.5)  characterized  by 

'0.1000  0.0490  0.0140  0.0105' 

0.0490  0.0700  0.0275  0.0130 

Q = 

0.0140  0.0275  0.1200  -0.0208 

0.0105  0.0130  -0.0208  0.0600 

_ - 

A standard  predictive  controller  is  designed  with  tuning  parameters  Ny  = 4, 

Nu  = 2 , and  A-  0 , resulting  in  the  following  controller  polynomials: 

R(z)  = z2  - 1.2920z  + 0.2920 

S(z)  = 1.7900 z2  - 0.9487z  + 0.1558 

T(z ) = -0.1578z4  - 0.0106z3  + 0.2905z2  + 0.8750z 

The  SPC  controller  stabilizes  the  nominal  system  because  the  characteristic  polynomial 

Aq(z)  = z4  - 0.302 lz3  has  all  its  roots  inside  the  unit  circle.  The  methodology  of 

Section  3.3.1  is  used  to  characterize  the  robustness  of  the  SPC  design,  yielding  a 

parametric  stability  margin,  a - 0.8833  < 1.  Hence,  this  design  is  not  robustly 

stabilizing,  and  it  can  be  concluded  that  there  are  plants  in  the  ellipsoidal  domain  for 

which  the  closed  loop  is  unstable. 
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Figure  3.2  illustrates  the  results  of  SPC  closed  loop  simulations  with  both  a 
nominal  and  perturbed  plant.  The  uncertain  coefficient  vector  considered  is 

<5p  = [0.0971  -0.0093  0.2470  -0.1 599] r 
which  corresponds  to  the  perturbed  plant 

B(z)  + Afl(z)  _ 1.2470z  - 0.4599 

A(z)  + A A(z)  z2  - 0.7029  + 0.1 507 

The  coefficient  vector  is  included  in  the  ellipsoidal  domain  because 
<5prQ<5p  = 0.9714  <1.  The  SPC  design  is  tested  via  servo  and  regulation  tests  by 
making  a setpoint  change  at  t = 75  and  introducing  a persistent  step  disturbance  d(t)  = 1 
at  t = 25.  The  SPC  stabilizes  the  loop  of  the  nominal  plant,  but  produces  an  unstable 
closed  loop  for  the  perturbed  plant,  causing  unbounded  input  and  output  signals. 
Therefore,  the  standard  p redictive  c ontroller  is  not  acceptable  from  the  r obust-stability 
viewpoint  and  it  must  be  redesigned. 


t 


Figure  3.2.  Ellipsoidal  uncertainty  example.  Closed  loop  responses  of  the 
SPC  design  considered  in  the  example  to  a set  point  change  and  to  an 
exogenous  disturbance  injection. 
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A robust  predictive  controller  (RPC)  is  designed  by  first  calculating  the  value  of 
the  robust  stability  margin,  a using  the  numerical  method  discussed  in  Section  3.1  over  a 
range  of  the  control-weighting  parameter,  namely  0 < X < 10 . The  prediction  and  control 
horizons  are  kept  at  the  same  values  used  by  the  SPC  design,  namely,  Ny  = 4 

and  Nu  = 2 . It  is  straightforward  to  verify  that  for  all  the  values  of  X considered,  the 

poles  of  Aq(z ) lie  strictly  in  the  unit  circle.  Therefore,  all  values  0<T<10  are 
nominally  stabilizing.  The  results  are  shown  in  Figure  3.3,  where  the  parametric  stability 
margin  initially  increases  sharply,  reaching  a maximum  value  of  a = 1.715  at  X = 0.6  . 
All  values  of  X in  Figure  3.3  that  result  in  a > 1 ensure  robust  stability.  Note  that 
Figure  3 .3  shows  that  all  designs  with  X < 0.0413  are  not  robustly  stabilizing.  U sing 
Figure  3.3,  the  weighting  parameter  is  selected  as  X = 0.6  for  maximum  stability 
robustness.  The  resulting  maximally  robust  RPC  design  is  then  tuned  using  X = 0.6 , 
Ny  - 4 , and  Nu  = 2 . 


Figure  3.3.  Ellipsoidal  uncertainty  example.  Effect  of  the  control- 
weighting parameter  X on  the  parametric  stability  margin  a . The  values  of 
a located  above  the  horizontal  dashed  line  correspond  to  robustly 
stabilizing  weighting-parameter  values. 
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Figure  3.4.  Ellipsoidal  uncertainty  example.  Closed  loop  responses  of  the 
RPC  design  considered  in  the  example  to  a set  point  change  and  to  an 
exogenous  disturbance  injection. 

The  RPC  design  stabilizes  the  nominal  plant  since  the  resulting  nominal 
characteristic  polynomial  Aq(z)  = z4  - 0.792 lz3  + 0.2724z2  - 0.0386z  is  Schur  given 
that  its  zeros  are  located  at  {0,  0.2380 ± 0.2556y  , 0.3161}.  The  RPC  design  is  also 
robustly  stabilizing  because  a - 1.715.  Figure  3.4  illustrates  the  RPC  simulation  results 
for  the  nominal  and  perturbed  plant  using  the  same  servo  and  regulation  tests  performed 
earlier.  Note  that  the  RPC  design  produces  adequate  output  responses,  effectively 
tracking  the  setpoint  and  rejecting  the  disturbance  for  both  the  nominal  and  the  perturbed 
plant. 

3.4.2  Hyperbox  Uncertainty 

Consider  the  nominal  plant  (3.28)  with  an  associated  hyperbox  uncertainty  of  the 


b = [0.2  0.2  0.2  0.2] T 


form  (3.8)  where 
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As  in  the  previous  example,  two  controller  designs,  denoted  as  SPC  and  RPC  are 
evaluated  for  robust  stability  and  performance. 

The  SPC  tuning  parameters  chosen  are  Ny  - 4 , Nu  - 2 , and  X = 0.1 . The  closed 

loop  under  SPC  control  is  nominally  stable  because  the  characteristic  polynomial 
Aq (z)  = z4  -0.4528z3  + 0.0831z2  -0.0116z  has  all  roots  inside  the  unit  circle.  Using 
the  linear-programming  method  discussed  in  section  3.3.2,  it’s  readily  determined  that  the 
robust  stability  of  the  SPC  design  has  a parametric  robust-stability  margin 
a - 0.5942  < 1 . Therefore,  the  SPC  design  is  not  robustly  stabilizing  with  respect  to  the 
hyperbox  uncertainty  considered. 

Given  that  a < 1 for  the  SPC  design,  it  is  anticipated  that  the  closed  loop  must  be 
unstable  for  some  plants  in  the  uncertainty  domain.  Consider  the  coefficient  vector 

£>p  = [0.1701  -0.1701  0.1701  — 0.1 701]  ^ in  the  hyperbox  domain,  which  results  in  the 
perturbed  plant 

B(z)  + Aff(z)  1.701z- 0.4701 
A(z)  + AA(z)  z2  - 0.6299  - 0.0101 

Figure  3.5  illustrates  that  the  SPC  controller  stabilizes  the  nominal  plant,  but  produces  an 
unstable  response  for  the  perturbed  plant.  A more  robust  controller  must  be  designed. 
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Figure  3.5.  Hyperbox  uncertainty  example.  Closed  loop  responses  of  the 
SPC  design  considered  in  the  example  to  a set  point  change  and  to  an 
exogenous  disturbance  injection. 

The  RPC  is  designed  taking  advantage  of  the  effect  of  A the  parametric  stability 
margin,  a , as  illustrated  in  Figure  3.6  where  the  horizons  are  fixed  at  N y = 4 and 

Nu  - 2 . It  is  straightforward  to  verify  that  all  values  of  A considered  (namely 
0<  X<  10)  produce  nominally  stabilizing  predictive  control  designs  ( i.e the  resulting 
Aq(z)  are  Schur).  The  parametric  stability  margin  increases  sharply  from  X = 0 to 
X = 2.7  where  it  realizes  its  maximal  value  a = 1.626 . Any  X that  results  in  a > 1 
ensures  robust  s tability  with  r espect  t o t he  given  u ncertainty  set.  The  choice  X - 2.7 
corresponds  to  the  RPC  design  with  maximal  robust  stability.  Note  that  all  weighting 
parameter  choices  X < 0.3825  fail  to  deliver  a robustly  stabilizing  design  because  a < 1 . 
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Figure  3.6.  Hyperbox  uncertainty  example.  Effect  of  the  control-weighting 
parameter  A on  the  parametric  stability  margin  a.  The  values  of  a located 
above  the  horizontal  dashed  line  correspond  to  robustly  stabilizing 
weighting-parameter  values. 
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Figure  3.7.  Hyperbox  uncertainty  example.  Closed  loop  responses  of  the 
RPC  design  considered  in  the  example  to  a set  point  change  and  to  an 
exogenous  disturbance  injection.. 

The  RPC  is  designed  using  the  tuning  parameters  Ny  = 4 , Nu  = 2 , and  A = 2.7 . 
Robust  stability  is  ensured  because  a = 1.626  > 1 . Nominal  stability  is  also  attained 
given  that  the  characteristic  polynomial  Aq(z)  = z 4 -1.1048z3  + 0.4588z2  -0.0673z  is 
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Schur  because  its  zeros  are  {0,  0.3849  ± 0.2297 j , 0.3349}.  Figure  3.7  illustrates  that  the 
RPC  stabilizes  the  closed  loop  featuring  either  the  nominal  or  the  perturbed  plant. 

3.5.  Conclusions 

This  chapter  presents  an  approach  to  analyze  robust  predictive  controllers  based 
on  a parametric  stability  margin  that  can  be  computed  numerically  using  well-known 
algorithms.  The  parametric  stability  margin  offers  a quantitative  measure  of  robust 
stability  that  destabilizes  the  system,  and  serves  as  a basis  for  the  efficient  robustification 
of  predictive  controller  designs.  Two  types  of  parametric  uncertainty  descriptions, 
namely  ellipsoidal  and  hyperbox,  are  treated  in  detail.  The  systematic  algorithms 
proposed  allow  the  user  to  determine  the  effect  of  a tuning  parameter  on  the  robust 
stability  of  the  system.  Although  the  examples  focused  on  the  effect  of  the  control- 
weighting parameter  X , an  analogous  procedure  can  be  used  to  analyze  the  effects  of  the 
control  horizons  Nu  and  Ny  . 


CHAPTER  4 

PREDICTIVE  CONTROL  OF  A WASTEWATER  TREATMENT  PLANT 

4.1  Introduction 

Government  regulations  require  wastewater  treatment  plants  to  disinfect  water 
before  it  is  discharged  to  the  environment.  This  can  be  accomplished  in  a variety  of 
methods,  but  the  most  common  process  is  to  add  chlorine  to  reduce  the  harmful 
compounds  (Metcalf  & Eddy,  1979).  Chlorine  reacts  in  this  process  and  regulations 
necessitate  a minimum  level  of  chlorine  concentration  and  minimum  residence  time  upon 
leaving  the  reactor.  This  ensures  that  enough  chlorine  is  present  to  provide  satisfactory 
cleansing  (Hammer  and  Hammer,  1996).  However,  high  chlorine  levels  are  harmful  for 
three  reasons:  (1)  chlorine  is  a significant  operating  cost,  (2)  unnecessary  addition  of 
chlorine  means  that  greater  volumes  must  be  stored  on-site,  and  (3)  high  chlorine 
concentrations  tend  to  form  trihalomethane  (THM)  contaminants  (Jolley  et  al.,  1990). 
Therefore  t he  c ontroller  d esign  c annot  a llow  h igh  c oncentrations  o f c hlorine  and  must 
maintain  the  outlet  chlorine  concentration  at  the  desired  set  point. 

The  objective  of  this  research  is  to  develop  an  effective  chlorine  control  scheme 
to  overcome  the  technological  challenges  associated  with  a wastewater  treatment  plant  . 
Plants  have  little  storage  space  for  water;  thus  they  must  process  water  as  it  arrives  and 
are  subject  to  large  cyclical  demands  for  water  (Hammer  and  Hammer,  1996).  This 
change  in  water  flow  causes  a time-varying  delay,  due  to  the  reactor’s  size  being  fixed, 
and  is  a considerable  challenge  for  process  control  of  wastewater  (Smith,  1971).  There 
are  other  lesser  disturbances  such  as  intense  solar  radiation  (causing  increases  in 
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volatilization  of  chlorine,  evaporation  of  water,  and  chemical  reaction  rates)  and  rainfall 
(causing  dilution)  that  can  affect  the  dynamics  of  outdoor  reactors. 

Several  feedback  control  strategies  have  been  developed  throughout  the  years  to 
provide  optimal  outlet  chlorine  concentration  levels  with  varying  results.  The  simplest 
method  is  applying  a proportional-integral  (PI)  controller  in  a cascade  feedback  form,  but 
this  usually  fails  for  systems  with  large  disturbances  and  time  delays  (Seborg  et.  al., 
1989).  The  addition  of  a Smith  Predictor,  requiring  an  accurate  model  of  the  system,  can 
reduce  the  effects  of  time  delays  if  the  delay  is  relatively  static  (Ogunnaike  and  Ray, 
1994).  These  two  methods  fail  to  control  a system  with  very  high  disturbances  or  time 
delays  as  are  present  in  wastewater  treatment  plants.  Meredith  (2003)  found  an  effective 
technique  by  modifying  the  reactor  with  a dynamic  weir  which  allows  the  controller  to 
manipulate  the  residence  time  of  the  water  in  the  reactor  by  changing  the  reactor  volume. 
This  modification  essentially  makes  the  time  delay  of  the  system  constant  and  has  been 
found  to  be  extremely  effective  in  controlling  chlorine  concentration.  However  this 
requires  a physical  modification  to  the  wastewater  treatment  plant. 

Many  advanced  control  strategies  require  accurate  models  of  the  process  in  order 
to  be  effective.  A schematic  of  the  general  design  of  wastewater  treatment  plants  is 
shown  below  in  Figure  4.1.  The  wastewater  enters  the  dosing  area  where  chlorine  is 
added  to  eliminate  ammonia  nitrogen  from  the  water  (Metcalf  & Eddy,  1979).  The 
ammonia  is  quickly  removed  in  this  basin  and  the  reactions  can  be  described  by  the 
Morris-Wei  model  (Morris  and  Wei,  1969).  However  the  Morris-Wei  model  introduces  a 
combination  of  algebraic  and  differential  equations  which  increase  the  computational 
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burden.  A simplification  was  developed  by  Snooeyink  and  Jenkins  (1980)  that  assumes  a 
first  order  model  and  greatly  simplifies  the  computation. 

Chlorine  Input 

► 

Wastewater 


0 L 

Figure  4. 1 . Schematic  of  wastewater  treatment  plant 

The  outflow  from  the  dosing  area  then  enters  the  contact  area  where  most  of  the 

disinfection  takes  place.  This  area  of  the  reactor  is  modeled  as  a series  of  continuous 
stirred  tank  reactors  (CSTR)  in  series.  For  more  information  on  the  model  please  see 
Meredith  (2003)  and  Bierck  et  al  (1996).  Meredith  found  that  the  optimal  number  of 
CSTRs  that  best  modeled  the  plug  flow  reactor  was  100.  The  model  that  was  used  in  this 
study  for  the  development  of  predictive  controllers  and  to  assess  the  performance  of  the 
different  designs  was  co-developed  by  Demir  (2001)  and  Meredith  (2003). 

This  chapter  reports  on  the  development  of  a disturbance  predictive  controller  that 
was  developed.  In  section  4.2  the  derivation  of  the  control  law  is  illustrated  with  all 
pertinent  equations.  Section  4.3  illustrates  the  difference  between  the  responses  of  a 
disturbance  predictive  controller  and  predictive  controller  employed  with  both  a linear 
and  nonlinear  model.  Finally,  4.4  details  some  future  directions  and  control  schemes  that 
may  be  useful. 

4.2  Derivation  of  Disturbance  Predictive  Controller 

The  derivation  of  the  disturbance  predictive  control  law  follows  a similar  form  as 
shown  by  Crisalle  et  al.  (1989).  Consider  the  nominal  plant  transfer  function 


Dosage  Area 

Contact  Area 

V 

c 

57 


y(t)  = - \)+<Mh.d(t) 

A(q~X)  D(q~X) 


(4.1) 


with  input  u(t  - 1) , output  y(t),  disturbance  d(t) , and  known  coprime  polynomials 
A(q _1)  and  B(q~X) . Multiply  both  sides  of  (4.1)  by  A(q~X)  and  D(q~ *)  to  yield 


D(q-l)A(q-l)y(t ) = D(q-X)B(q-X)u(t  - 1)  + A{q-X)C{q-X)d{t)  (4.2) 


and  multiply  (4.2)  by  Et{q  X)Aql 

Ei(q~{)D(q~X)A(q~l)Ay(t  + i)  = EAq~X)D(q~')B(q~l)Au(t  + i-\)  (43) 

+ Ei(q-X)A(q-')C(q-l)Ad(t  + i) 

where  q 1 is  the  forward  shift  operator  and  A = 1 -q~X . Define 

Ad(q-X)  = D(q-X)A(q-X) 
and  consider  the  following  Diophantine  equation 

\ = Ei(q-')Ad(q-')t±  + q-iFi{q-1) 

The  Diophantine  equation  is  used  to  separate  the  future  values  of  y(t  + i ) from  the  past 
values  stored  in  y(t) . Rearranging  the  above  equation  yields 


Et(q-')A<I(q-')^  = \-q-iFl(q-')  (4.4) 

Equation  (4.4)  is  substituted  into  (4.3)  to  separate  y(t)  and  y{t  + /)  resulting  in 

y(t  + i)-Fi(q-l)y(t)  = Ei(q-X)D(q-')B(q-X)Au(t  + i-l)  + Ei(q-l)A(q-')C(q-l)Ad(t  + i)(4.5) 
Next,  the  past  and  future  values  of  the  input  must  be  similarly  separated  using  the 
Diophantine  equation  by  introducing  the  following  polynomials 

ri  (d~l)  = Yi,0  + Yi,  l?'1  +'  • ■+Yi,n^~Tlb  +1 

Qiq  x)  = g0+g\q  l+"-+gi-i<i  *+1 
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and  requiring  V t(q  *)  and  Gt(q  !)  to  satisfy 

Ei(q-l)Bd(q-l)  = Gi(q-l)  + q-iri(q-1)  (4.6) 

where 

Bd(q-l)  = D(q-l)B(q-1) 

To  separate  the  input  values  substitute  (4.6)  into  (4.5) 

y(t  + i)  = G'iq-^Auit  +i-l)  + ri(q-')Au(t  -l)  + Eiiq-^Aiq-'Mq-^Adit  +i)  + Giq-'Mt) 
and  make  the  following  definitions  for  yOL ( t+i ) and  yd (t  + i ) 

y°L(t  +i)  = Tiiq-^Auit  - 1)  + Ft(q~l)  (4.7) 

yd(t+i)  = Ei(q-l)A(q-1)C(q-1)Ad(t+i)  (4.8) 

The  equations  for  yOL{t  + i)  and  yd(t  + i ) are  utilized  later  in  the  development  of  the 
control  law. 

The  disturbance  predictive  control  law  is  derived  from  the  following  cost  function 

Nr  2 N 

J(u)=YJ[r(t  + i)-y(t  + i\t)\  +/t£[AM(* + 0]  (4.9) 

j=1  i=0 

where  r(t  + i ) are  the  future  values  of  the  set  point,  y(t+i\t)  are  the  predicted  future 
values  of  the  system  evaluated  at  time  t,  and  A u{t  + i)  are  the  future  values  of  the  inputs. 
Rewriting  (4.9)  in  vector-matrix  form  yields 

J(u)  = (y-r)r(y-r)  + Auru  (4.10) 

Define  the  output  of  the  system  as 

y = GAu  + yOL  +yd 


and  substitute  into  (4.10) 
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T 

J(u ) = (GAu  + yOL  + yd  - r)  (GAu  + yOL  + yd  - r)  + Au  Tu 


Define  the  initial  input  as 


uo  = r-y 


OL 


and  simply  (4.1 1)  to  yield 

T 

J(u ) = (uq  - GAu  - yd ) (uq  — GAu  - yd  j + Jlu  Tu 

Take  the  derivative  of  equation  (4.12)  to  produce 

^^  = -UoG  + u7’G7’G  + ydrG  + Aur  =0 
8u 

Solving  (4.13)  for  u 


u = 


G G + 2I 


r-yOL  -yd 


and  rewriting  yields 


A u{t)  = k 


w(t)-yOL{t  + i)-yd(t  + i ) 


where 


w(t  + i)  - qi r(t) 


k =[10---0]K 


and 


K = Gl 

where  G \ is  the  truncated  dynamic  matrix. 

Next,  put  equation  (4.14)  in  summation  form 

Nu  N,  Nt  Nt 

Z Am(? + 0 = Z kM{ + 0 - Z kiyOL  (* + 0 - Z ktyd  0 + 0 

i=0  i'=l  i=l  i'=l 


(4.11) 


(4.12) 


(4.13) 


(4.14) 


(4.15) 
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and  substitute  (4.7)  and  (4.8)  into  (4.15)  to  yield  the  control  law 

N 

i+r'zw?-1) 

i=l 

-'ZkiEi(q-l)A(q-l)C(q-1)Ad{t  + i) 

i=l 

The  control  law  described  by  equation  (4.16)  can  be  expressed  as  R(z) , S(z ) , T(z) , and 
V (z)  polynomials  defined  as 


N 


N 


a u(t)  = 2>^V(0-  !)t(0 


i=i 


1= 1 


R(z)  = z" 


\+z~'Y.kiVi(z-') 

J = 1 


N 


-1> 


A 


S(z) = zn 


’N 

IW‘) 

1=1 


N 

T(z)=ikizi 

i=i 


JV 


n*)  = EW*-1)^-1)^"1)* 


j=i 


and  is  illustrated  in  graphical  form  below  in  Figure  4.2. 
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Disturbance 


Figure  4.2.  Structure  of  a disturbance  predictive  controller  characterized  by 
the  operator  polynomials  R(z ) , S(z) , T(z) , and  V(z) . 


4.3  Simulation  of  Wastewater  Treatment  Plants 

In  this  section,  the  predictive  control  law  developed  by  Crisalle  et  al.  (1989)  is 
compared  to  the  disturbance  predictive  controller  developed  in  section  4.2.  Each 
controller  is  connected  to  the  model  developed  by  Meredith  (2003)  and  the  closed  loop 
response  obtained.  In  section  4.3.1,  a linear,  first  order  plus  deadtime  model  is  used  as  an 
approximation  of  the  highly  nonlinear  wastewater  treatment  plant  and  utilized  as  an 
initial  test  system  for  the  comparison  of  the  two  controllers  is  illustrated.  The  controller 
parameters  discovered  during  the  linear  model  simulations  are  used  to  control  the 
nonlinear  model  of  the  wastewater  treatment  plant  discussed  in  section  4.3.2 
4.3.1  Closed-Loop  Response  of  Linear  Model 

Figure  4.3  shows  the  performance  of  the  predictive  control  system  developed  by 
Crisalle  et  al.  (1989)  applied  to  the  linear  plant  model  when  the  flowrate  is  modeled  as 
sinusoidal  wave,  which  is  chosen  to  approximate  a smooth,  cyclical  water  consumption 

profile  of  a typical  day.  The  sine  wave  begins  at  time  6 hr,  its  period  is  24  hr,  and  the 

3 3 

amplitude  is  three-fourths  of  the  nominal  flowrate  (0.4  m /s),  or  0.3  m / s.  The  result  is  a 
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disturbance  flowrate  that  varies  from  0.1  m /s  to  0.7  m /s.  These  values  are  chosen 

because  they  represent  a 7:1  ratio  of  maximum-to-minimum  flowrate. 

The  linear  plant  is  controlled  with  unacceptable  performance,  experiencing 
fluctuations  in  the  output  chlorine  concentration.  The  sinusoidal  flow  of  the  amplitude 
described  above  severely  degrades  the  quality  of  the  controller  performance  for 
simulations  of  the  linear  plant.  Figure  4.3  shows  that  the  deviation  output  concentration 
values  are  as  high  as  1.27  ppm  and  as  low  as  -1.45  ppm  which  would  correlate  to  3.27 
ppm  and  0.55  ppm  if  the  set  point  was  2 ppm.  The  lower  value,  0.55  ppm  of  chlorine,  is 
unacceptable  because  it  is  below  the  level  of  1 ppm  established  by  the  EPA. 
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Figure  4.3.  Predictive  control  of  linear  wastewater  treatment  plant  to  a 
sinusoidal  disturbance 
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The  disturbance  predictive  controller  developed  in  section  4.2,  however,  does  an 
excellent  job  of  maintaining  the  outlet  chlorine  concentration  as  illustrated  in  Figure  4.4. 
The  sinusoidal  flow  of  the  amplitude  described  above  does  not  degrade  the  quality  of  the 
controller  performance  for  simulations  of  the  linear  plant.  Figure  4.4  shows  that  the 
deviation  output  concentration  values  are  only  as  high  as  0.29  ppm  and  as  low  as  -0.26 
ppm  which  would  correlate  to  2.29  ppm  and  1.74  ppm  if  the  set  point  was  2 ppm.  Both 
of  these  values  are  well  within  the  regulations  established  by  the  EPA  demonstrating  the 
improvement  i n c ontroller  p erformance  i f t he  d isturbance  p redictive  control  s cheme  i s 
utilized. 
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Figure  4 .4.  D isturbance  p redictive  c ontrol  o f 1 inear  w astewater  t reatment 
plant  to  a sinusoidal  disturbance 
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4.3.2  Closed-Loop  Response  of  Nonlinear  Model 

The  physical  wastewater  treatment  plant  is  a highly  nonlinear  process  with  very 
different  dynamics  than  those  seen  in  the  linear  model.  In  Figure  4.5  a disturbance 
predictive  controller  is  connected  to  the  nonlinear  model  of  the  wastewater  treatment 
plant.  The  sinusoidal  disturbance  is  the  same  as  used  in  section  4.3.1  and  results  in  an 
outlet  value  of  chlorine  concentration  of  0.3  ppm  and  5.9  ppm.  These  upper  and  lower 
values  are  not  allowable  by  the  EPA  rendering  the  disturbance  predictive  controller 
unacceptable  for  use  in  the  highly  nonlinear  wastewater  treatment  plants. 


0 10  20  30  40  50  60 


Figure  4.5.  Closed-Loop  response  of  nonlinear  m odel  with  a disturbance 
predictive  controller  subjected  to  a varying  flow 

4.4  Conclusions  and  Future  Work 


The  disturbance  predictive  controller  developed  is  similar  in  structure  and  nature 
to  previous  predictive  controllers.  However  it  needs  an  accurate  model  of  the  system  in 
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question  as  well  as  the  disturbance  in  order  to  be  effective.  The  controller  is  linear  in 
form  and  does  an  excellent  job  if  the  system  to  be  controlled  is  linear  as  illustrated  in 
section  4.3.1.  However  if  the  system  is  highly  nonlinear,  such  as  a wastewater  treatment 
plant,  with  very  large  disturbances,  the  disturbance  predictive  controller  fails.  One 
solution  is  to  design  an  adaptive  disturbance  predictive  controller  similar  in  form  to  that 
described  by  Landau  (1998).  However,  this  would  be  very  numerically  intensive  as  the 
control  law  polynomials  need  to  be  updated  often.  A more  practical  and  effective 
approach  would  be  to  develop  a nonlinear  predictive  controller.  This  would  allow  the 
controller  to  maintain  a strict  chlorine  concentration  despite  the  large  disturbances. 


CHAPTER  5 

AN  EXCESS  VOLUME  MODEL  TO  DETERMINE  CONVERSION  IN  AN 
EMULSION  POLYMERIZATION  REACTOR 

5.1  Introduction 

Emulsion  polymerization  became  commercialized  in  the  1930s  and  has  grown  to 
become  the  major  process  for  the  production  of  synthetic  polymers  (Blackley,  1975). 
Much  of  this  development  was  generated  by  the  loss  of  natural  rubber  during  World  War 
II.  Postwar,  the  wide  spread  usage  of  water  based  paints  expanded  emulsion 
polymerization  processes  to  include  acrylic  resins  and  synthetic  elastomers.  Presently, 
many  different  types  of  polymers  are  produced  by  emulsion  polymerization  primarily 
through  batch  reactors.  However,  continuous  stirred  tank  reactors  are  more  commonly 
utilized  due  to  economic  and  process  advantages  (Gilbert,  1995). 

There  is  a lack  of  data  on  emulsion  polymerization  reactors  available  in  literature. 
This  slow  development  of  polymerization  reactors  is  due  to  the  difficulty  in  developing 
sensors  that  can  monitor  important  polymer  properties  on-line  (Schork  and  Ray,  1980; 
Schork  and  Ray,  1987). 

Density  is  measured  on-line  through  the  use  of  a densitometer,  a vibrating  U- 
shaped  tube  which  continuously  has  liquid  flowing  through  it.  The  first  important  work 
concerning  on-line  densitometry  was  performed  by  Canegallo  et  al.  (1993)  where  thermal 
instability  and  coalescence  of  monomer  droplets  was  analyzed  and  a phase  separator 
introduced  before  the  sample  stream  enters  the  densitometer.  This  phase  separator 
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improves  density  measurements  by  removing  the  gas  from  the  liquid  sample  to  be 
analyzed. 

Schork  and  Ray  (1981,  1983)  were  the  first  to  measure  conversion  on-line  by 
developing  an  equation  that  related  density  to  conversion  with  an  emphasis  on  emulsion 
polymerization.  They  assumed  that  in  an  emulsion  polymerization  reaction  all  volumes 
of  the  components  were  additive  and  verified  their  experimental  results  with  gravimetric 
analysis.  Ponnuswamy  (1986)  improved  the  equation  by  comparing  densitometry  and  GC 
analysis  of  data  and  was  able  to  greatly  improve  the  estimate  of  the  conversion. 

A large  effort  is  currently  underway  to  development  on-line  sensors  to  monitor  the 
emulsion  polymerization  process.  Guinot  et  al.  (2000)  performed  a comparative  study 
and  analyzed  the  strengths  and  weaknesses  of  many  sensors  such  as  IR,  NER,  and  polymer 
reaction  calorimetry.  Vara  (2000)  improved  the  densitometer  measurements  and  analyzed 
such  issues  as  thermal  instability  for  a wide  variety  of  polymerization  processes.  This 
chapter  takes  the  work  begun  at  the  University  of  South  Florida  by  Vara  and  improves  the 
conversion  calculations  through  the  use  of  an  excess  volume  model.  Numerical  issues  in 
root  finding  are  addressed  and  several  examples  are  presented. 

5.2  Motivation 

A realistic  model  for  describing  the  specific  volume  of  a mixture  of  polymer  and 
its  monomer  in  an  emulsion  polymerization  reactor  is  of  the  form 

v W = vo%  + (V100%  - vo%  )x  + V-  (*)  (5.1) 

(x)  = v(x)  - (v0%  + (v100%  - v0%  )*) 

where  x is  the  degree  of  conversion  of  monomer,  v(x)  is  the  specific  volume  of  the 
mixture  at  the  given  conversion  level  x,  v0%  is  the  specific  volume  of  the  monomer  at  0% 
conversion,  v/oo%  is  the  specific  volume  of  polymer  at  100%  conversion,  and  vex(x)  is  the 
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excess  volume  of  mixing.  Equation  (5.1)  can  be  interpolated  as  the  thermodynamic 
equation  for  the  volume  of  a mixture,  and  is  based  on  the  Lewis-Randall  rule  to  select  the 
standard  states  as  vo%  and  V]0o%. 

The  specific  volume  of  an  emulsion  can  be  determined  by  measuring  the  density  p 
of  a sample  and  then  invoking  the  relationship 

v(*)  = - 

P 


Given  that  the  variable  v(x)  can  be  measured,  in  principle,  it  is  possible  to  solve  (5.1)  for 
the  remaining  unknown  variable,  namely,  the  degree  of  conversion  x.  Experimental  work 
carried  out  by  Vara  assumed  that  the  excess  volume  of  mixing  is  negligible  (i.e., 
vex(x)=0),  and  succeeded  in  obtaining  experimental  values  of  degree  of  conversion 
through  the  expression 


x = 


v(*)-v, 


0% 


V100%  vo% 


(5.2) 


Vara  compared  the  results  obtained  in  this  fashion  with  the  trustworthy  degree-of- 
conversion  values  obtained  via  gravimetric  analysis,  and  found  that  the  approach  (5.2) 
underestimates  the  degree  of  conversion.  The  conclusion  that  can  be  inferred  from  these 
observations  is  that  it  is  necessary  to  take  into  account  the  excess  volume  of  mixing  term 
in  (5.1)  for  improved  accuracy  in  the  estimation  of  the  degree  of  conversion. 

An  example  of  the  motivation  behind  the  excess  volume  modeling  is  illustrated  in 
Figure  5.1.  The  linear  relationship  between  conversion  and  specific  volume,  as  shown  by 
equation  (5.2),  is  plotted  and  called  the  ideal  model  (v“l(x)).  The  actual  specific  volume  is 
shown  b y t he  c urve  v (x)  and,  in  this  case,  lies  above  the  ideal  model.  Therefore,  if  a 
specific  volume  measurement  is  made,  the  ideal  model  will  result  in  a much  smaller 
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conversion  value  than  is  present  in  the  reactor.  Consequently,  an  excess  volume  model 
will  significantly  improve  conversion  calculations. 


vioo% 


Figure  5.1.  Motivation  for  use  of  an  excess  volume  model 
In  this  chapter,  a one-parameter  model  and  two-parameter  model  for  the  excess 

volume  of  mixing  are  proposed.  The  choice  of  this  model  depends  on  the  shape  of  the 

excess  volume.  For  example,  if  the  excess  volume  is  symmetric,  a one-parameter  model 

may  be  used. 


5,3  One-Parameter  Excess  Volume  Model 
5.3.1  Mathematical  Formulation 

For  a symmetric  excess  volume,  the  following  one-parameter  model  is  proposed: 

v“(x)  = yl,x(l-x)  (5.3) 

where  the  parameter  Ai  must  be  determined  experimentally.  The  parameter  A/  can  be 
defined  as  the  excess  volume  coefficient.  The  one-parameter  model  assumes  that  the 
dependence  of  vex(x)  on  x is  symmetric  in  the  conversion  range  [0,  1];  however,  this 
assumption  must  be  verified  experimentally.  If  through  further  experimental  studies  it  is 
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determined  that  excess  volume  curve  is  asymmetric,  then  a two-parameter  model  can  be 
considered. 

From  (5.3)  it  follows  that  the  maximum  absolute  value  of  the  excess  volume  is 
realized  at  50%  conversion  (i.e.,  x=0.5).  Defining  v^%:=  vet(0.5)  , it  follows  that 

vs“,,  =40.5(1-0.5)  = A 

Or,  equivalently, 

A =4v5£o% 

Therefore,  the  modeling  parameter  A / can  be  interpreted  as  the  excess  volume  of  mixing 
at  50%  conversion  multiplied  by  a factor  of  4. 

Note  that  vex(x)  (and  therefore  the  parameter  A/)  must  be  a function  of  the 
molecular  weight  distribution  of  the  polymer.  Analogously,  the  variable  v/oo%  also 
depends  on  the  molecular  weight  distribution.  Throughout  these  developments  it  is 
assumed  that  both  A/  and  vioo%  are  weak  functions  of  the  molecular  weight  distribution, 
as  is  implicitly  assumed  by  equation  (5.2). 

Substituting  the  model  (5.3)  into  the  specific  volume  expression  (5.1)  yields  the 
corrected  specific  volume  equation 

v(x)  = v0%  + (v100o/o  — v0o/o  )x  + Aix(\  — x)  (5 .4) 

or 

N-X-)  = V0%  + (V100%  — V0%  + At)x  — AjX 
5.3.2  Determining  the  Excess  Volume  Parameter 

The  excess  volume  parameter  Aj  can  be  obtained  from  an  experimental  value  of 
conversion  x,  determined  from  gravimetric  analysis  and  from  the  corresponding  value  of 
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specific  volume  v(xi)=l/p(xI)  obtained  using  a densitometer  measurement.  Substituting 
both  of  these  experimental  values  in  (5.4)  and  rearranging  terms  gives  the  expression 


which  can  be  readily  solved  for  the  unknown  to  yield 

^ vo%  ~ v(xi ) + (vioo%  ~ vo%  )-xi  , (5  6) 

The  experimental  conversion  value  xy  must  be  selected  away  from  the  extreme  values  of 
x=0  andx=l  where  excess  volume  model  (5.3)  vanishes,  rendering  equation  (5.6)  invalid. 

A better  method  for  obtaining  an  experimental  value  foryl  y is  by  least-squares 
regression  using  a finite  number  of  experimental  conversion  data  x,  and  a corresponding 
set  of  specific  volume  data  v(xd,  /= 1,  2,  ...,  N.  As  discussed  earlier,  the  experimental 
points  should  be  chosen  far  away  from  the  extreme  values  of  x,=0  and  xt=\.  The  data  are 
then  used  to  pose  the  following  set  of  A equations  patterned  after  equation  (5.5): 


Finally,  the  system  of  equations  is  readily  solved  in  the  least-squares  sense  to  give  the 
parameter  estimate 


y4,X,(l  X,)  — Vo%  v(-*l)  + (V100%  V0%)*1 


(5.5) 


— Axx  j (1  Xj ) — v0%  v(x, ) + (v100o/(i  v0%  )Xj 

— Axx2  (1  — x2 ) = v0o/o  — v(x2 ) + (v100%  — v0o/o  )x- 


100% 


—AtxN(l  xN)  — v Oo/o  v(xN ) + (Vi00o/0  v0%)xw 


A,=- 


(5.7) 


Clearly,  equation  (5.7)  reduces  to  (5.6)  when  N=  1. 
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5.3.3  Solving  the  One-Parameter  Quadratic  Equation 

The  specific  volume  equation  (5.4)  can  be  solved  for  the  degree  of  conversion  x 
when  an  experimental  value  of  specific  volume  is  available.  Rearranging  terms  in  (5.4) 
yields  the  second-order  polynomial  equation 

Ax2  ~ (vi<x>%  - vo%  + A )*  + (Ax)  - v0% ) = o 

with  roots 


^.(1)  _ V100%  Vq%  + A y^V/o  V100%  Ay  4 A (V(x)  V0%) 
2 A,  2 A, 


(2)  _ V100%  V0%  + A V (V0%  V100%  A)  (v(x)  vm ) 

2 A,  2A, 


or 


r(»  - 


X2)  _ 


1 

— + 
2 


V1 00%  vo% 


2 A 


i y 


■ + 


V0%  V100% 


2 A. 


v(x)-v, 


0% 


1 / 


A 


( . _ A 

+ V10Q%  _ V0% 


2 A 


+ 


i 7 


J_  , V0%  V100% 


2 A. 


v(x)-v, 


0% 


(5.8) 


(5.9) 


i 7 


Therefore,  the  degree  of  conversion  x will  be  given  by  either  x=x<l>  or  x=x(2). 

Further  studies  are  needed  to  determine  how  one  should  select  one  of  the  two 
available  roots  as  a solution  to  the  conversion  problem.  One  case  where  the  problem  is 
easily  resolved  is  when  one  of  the  roots  lies  inside  the  interval  [0,  1]  and  the  other  root 
lies  outside.  In  this  case  the  selection  is  straightforward  because  the  solution  must  satisfy 
0 < x < 1 ; hence,  one  simply  chooses  the  root  that  lies  inside  the  [0,  1]  interval. 

A second  case  that  is  more  complicated  arises  should  both  roots  lie  inside  the 
interval  [0,  1].  Consider  for  example  the  case  where  a density  value  p=pioo%  is  measured 


at  an  instant  when  the  reactor  has  reached  100%  conversion.  It  then  follows  that 
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v(x)-VpIOoo%=vioo%;  substituting  this  quantity  into  (5.8)  and  (5.9)  and  carrying  out  simple 
manipulation  yields 


or 


x(1)  = 


x(2)  = 


1 

— + 
2 


V100%  vo% 


2 A 


i 


1 

— + 
2 


V0%  V100% 


2 A, 


V 

) 


, Vi  00%  vo%  ^ 


2 A 


+ . 


^ 1 _ V 

_|_  V0%  V100% 


1 J 


V 


2 A 


i J 


(5.10) 


(5.11) 


jc(1)  = 


}_  V100%  vo% 


2 A 


i / 


1 ! vo%  vioo% 


2/4, 


V0%  V100% 


1 J 


x(2)- 


J_  + ^100%  vo% 


2/4, 


+ 


i 2 


1 V1 00%  vo% 


2/4. 


= 1 


i / 


If  0 <(yo%  -vioo%)l  Aj<  1,  then  it  would  be  difficult  to  discern  between  the  correct  answer 
x=x<2>=\  and  the  alternative  answer  x=x<lj=  (v0%  -vioo%)l  At  because  both  roots  lie  in  the 
interval  [0,  1],  Further  research  must  be  carried  out  to  determine  if  the  condition 
0<(v0%  -vjoo%)/  Aj<  1 is  physically  meaningful  (for  example,  this  case  would  not  be 
relevant  if  physical  arguments  c an  b e u sed  t o s how  t hat  ( v0%  - vioo%)ly<  0 , e ffectively 
eliminating  the  root  x(1)).  A more  comprehensive  study  of  the  roots  in  this  case  will  be 
necessary,  particularly  taking  into  consideration  the  sign  of^y. 

A third  case  of  relevance  is  a situation  where  both  roots  x(\)  and  x(2)  lie  outside 
the  interval  [0,  1],  Such  roots  are  meaningless,  but  such  a situation  could  conceivably 
arise  if  the  densitometer  produces  an  incorrect  measurement  located  outside  the  range 


[ pi oo%,  Po%\  ■ Although  this  outcome  is  hypothetical,  it  should  be  ruled  out  as  a possibility 
to  e nsure  t he  r obust  p erformance  o f t he  d egree-of-conversion  p rediction.  P erhaps  t his 
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problem  can  be  dealt  with  easily  by  introducing  upper  and  lower  saturation  values  on  the 
acceptable  output  of  the  densitometer. 

A final  concern  that  must  be  taken  into  account  is  that  for  systems  where  A y is 
small  (i.e.,  the  excess  volume  is  negligible)  the  formula  (5.2)  must  be  used  in  lieu  of 
formulas  (5.10)— (5.1 1).  Note  that  the  equations  (5.1 0)-(5. 1 1)  fail  in  this  case  due  to  the 
appearance  of  the  parameter  A[  in  the  denominator,  causing  unacceptable  division-by- 
zero  errors  when  the  parameter  is  either  zero  or  very  small. 

5.4  Redlich-Kister  Two-Parameter  Model 
5.4.1  Mathematical  Formulation 

If  the  experimental  excess  volume  is  determined  to  be  asymmetric,  a two- 
parameter  model  must  be  implemented.  The  following  model  is  proposed: 

v“  (x)  = Ax(\  - jc)  + Bx(  1 - x)(2x  - 1)  (5.12) 

where  the  parameters  A and  B must  be  determined  experimentally.  This  model  does  not 
assume  that  vex(x)  must  be  symmetric  in  the  conversion  range  [0,  1]. 

The  two-parameter  model  adds  another  level  of  complexity  in  the  determination 
of  the  model  coefficients  and  conversion  values.  There  are  no  useful  analytical 
expressions  for  determining  A and  B as  there  were  with  the  one-parameter  model. 
However,  the  same  assumptions  (that  A,  B , and  v/oo%  are  weak  functions  of  the  molecular 
weight  distribution)  are  used  in  the  Redlich-Kister  model  as  well. 

Substituting  the  two-parameter  model  (5.12)  into  the  specific  volume  expression 
(5.1)  yields  the  corrected  specific  volume  equation 


v(x)  - v,,*  + (v100o/o  v0%)x  + Ax(\  x)  + Bx{\  x)(2x  1) 


(5.13) 
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5.4.2  Determining  the  Excess  Volume  Parameters 

Using  the  Redlich-Kister  model  adds  a degree  of  difficulty  to  the  analysis.  The 
coefficients  A and  B have  no  simple  analytical  solutions.  Instead,  least-squares 
regression  must  be  used  and  be  of  the  form  y-  Ax  with  a solution  x = (A'  A)~]  Ary 
where 


y- 


Lw 


A = 


*,(!-*,)  x1(l-x,)(2x1-l) 


x = 


_xn(\-xn)  xn(l-xn)(2xn-l)_ 

Therefore  if  B= 0 then  the  Redlich-Kister  model  reduces  to  the  one-parameter  model. 

5.4.3  Solving  the  Redlich-Kister  Cubic  Equation 

The  specific  volume  equation  (5.13)  can  be  solved  for  the  degree  of  conversion  x 
when  the  specific  volume  is  available.  A third-order  polynomial  results  when  the  terms 
in  (5.13)  are  rearranged 


x3  + 


35  A ) _2  | r A + v100o/o  Vq%  5 ^ | 


2 B J 


= 0 


2B  J V 2 B 

The  three  roots  for  this  expression  can  be  numerically  solved  to  yield  the  degree  of 
conversion. 

Further  study  is  needed  to  determine  how  to  select  one  of  the  three  roots  as  a 
solution  to  the  problem.  The  easiest  case  is  when  one  roots  lies  in  the  interval  [0,  1]  and 
the  other  two  lie  outside  the  interval.  In  this  case,  the  root  in  the  interval  [0,  1]  is  chosen. 

A second  case  that  is  more  complicated  occurs  when  one  of  the  roots  is  real  and 


two  are  complex  conjugates  of  each  other.  This  occurs  when  D> 0 where 
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D _ 1 (6^v10Qo/o  - 6Bvqq/o  + B2  + A2)3  | 


46656 
(3 B - A) 


B 1 


i2 


(5.14) 


24  B 


J~(A  + v100%  _ v0%  " B)  + — ^ + ~ 


2165- 


Therefore,  there  are  two  manners  in  which  D> 0.  First,  because  of  the  squared  term  in 
equation  (5.14),  D will  be  greater  than  zero  if 


{3B~A)(  u v v m i v°%  ~ vW  i (3jS  “ A) 

( A + v100%  - v0%  -B)+ 7^ + 3 

21653 


245^ 


> 


4 B 


1 (65v10o%  - 65v0o/o  + B2  + A2  )3 


46656 


5C 


which  will  guarantee  that  there  is  only  one  root.  An  alternative  condition  involves 
evaluating  the  cubed  term  in  equation  (5.14).  In  order  for  D> 0,  the  cubic  term  in  (5.14) 
must  be  negative.  Thus,  the  following  condition  must  hold 


6Bvm%  - 6Bv0%  +B~  + A"  <0 


Rearranging  terms  and  solving  for  vo%  -vioo%  results  in  the  expression 

A2  + B2 


6 B 


< V0%  V100% 


(5.15) 


From  the  physical  analysis  of  the  emulsion  system,  it  can  be  shown  that  the  quantity  vo%- 
vioo%  is  positive.  Therefore  if  the  A and  B coefficients  adhere  to  the  preceding 
expressions,  solving  for  the  degree  of  conversion  will  result  in  only  one  real  root. 

5.5  Example  - Polymerization  of  Styrene 


5.5.1  Development  of  Models 

An  experimental  emulsion  polymerization  apparatus  at  the  University  of  South  Florida 
was  used  to  o btain  g ravimetric  a nd  d ensity  d ata  v ersus  t ime  f or  t he  p olymerization  o f 
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styrene.  The  gravimetric  data  are  used  as  the  “correct”  value  for  the  conversion  in  the 
reactor.  Using  the  density  data  obtained  versus  time  and  equation  (5.2),  the  uncorrected 
value  of  conversion  was  calculated  and  compared  to  gravimetric  results  in  Figure  5.2. 


0 500  1000  1500  2000  2500 

t [s] 


Figure  5.2.  Comparison  of  gravimetric  and  densitometer  conversion  versus 
time 

As  illustrated  in  Figure  5.2,  the  gravimetric  and  densitometer  conversions  diverge  by  as 
much  as  0.15  from  0.4  to  0.8  conversion.  An  excess  volume  model  must  be  used. 

The  experimental  excess  volume  is  calculated  by  solving  equation  (5.1)  for  vex(x). 
A one-parameter  model  is  then  necessary  to  improve  the  densitometer  conversion 
calculations.  A least-squares  regression  is  performed  to  obtain  the  excess  volume 
parameter  A/,  as  shown  in  section  5.2.2.  The  experimental  excess  volume  is  then 
compared  to  the  one-parameter  model  for  vex(x)  in  Figure  5.3.  The  one-parameter  model 
is  symmetric  and  reaches  a maximum  at  0.5  conversion  while  the  experimental  vex  is  not 
symmetric  and  reaches  a maximum  around  0.75  conversion.  A two-parameter  model 
would  yield  a better  result  due  to  the  asymmetric  excess  volume. 
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In  order  to  obtain  the  parameters  A and  B,  a least-squares  regression  is  performed 
as  illustrated  in  Section  5.3.2.  The  Redlich-Kister  model  is  shown  in  Figure  5.3  and  is 
compared  to  the  experimental  excess  volume  and  the  one-parameter  model.  The  two- 
parameter  model  is  a better  fit  for  the  asymmetric  experimental  excess  volume.  However, 
the  model  still  predicts  a higher  excess  volume  between  0.3  and  0.6  and  underestimates 
the  excess  volume  between  0 and  0.3  as  well  as  0.6  to  0.9. 


-0.2  0 0.2  0.4  0.6  0.8  1 1.2 

X 


Figure  5.3.  Experimental  vex  and  vex  one-parameter  model  versus  conversion 

5.5.2  Conversion  Calculations 

In  Figure  5.4,  the  conversion  values  are  calculated  for  both  the  one-  and  two- 
parameter  models  and  plotted  versus  time.  Both  excess  volume  models  improve  the 
conversion  calculations;  however,  there  are  still  questions  associated  with  each. 

There  are  still  some  issues  concerning  the  one-parameter  model.  First,  at 
conversions  from  0.2  to  0.6,  the  c onversion  c alculated  for  the  o ne-parameter  m odel  i s 
actually  greater  than  for  the  gravimetric  conversion.  This  is  because  the  vex(x)  model 
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yields  a higher  excess  volume  than  is  shown  through  experimentation.  In  addition,  the 
model  predicts  less  conversion  between  0.6  and  0.8  than  is  reported  through  gravimetric 
determination.  A two-parameter  model  may  be  able  to  correct  these  inaccuracies. 

The  Redlich-Kister  model  is  then  used  as  the  excess  volume  model  and  the 
resulting  cubic  equation  is  solved  for  the  conversion.  Two  of  the  roots  must  be 
eliminated  in  order  to  solve  for  the  proper  value  of  conversion.  For  the  polymerization  of 
styrene,  equation  (5.15)  holds  true,  making  two  of  the  roots  imaginary.  The  resulting  real 
root  is  the  value  for  the  conversion  and  is  compared  to  the  one-parameter  model  and 
gravimetric  conversion  in  Figure  5.4. 


Figure  5.4.  Conversion  curves  for  different  excess  volume  models 
compared  to  gravimetric  and  densitometer  conversion 

As  illustrated  above,  the  two-parameter  model  provides  a more  accurate  value  for 


conversion  than  the  one-parameter  model  or  the  densitometer.  However,  at  conversions 
around  0.4  to  0.5,  both  the  one-  and  two-parameter  model  report  a conversion  higher  than 
the  actual  conversion  determined  by  gravimetric  methods  for  reasons  discussed 
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previously.  Possible  solutions  are  to  add  another  parameter  which  would  lead  to  a fourth- 
order  equation  in  terms  of  conversion  and  may  not  be  practical  to  solve. 

5.6  Copolymerization  of  Styrene  and  Butyl-Methacrylate 

5.6.1  Development  of  Models 

The  excess  volume  model  is  determined  by  taking  the  difference  between  the 
gravimetric  and  densitometer  conversion  at  different  conversions.  The  experimental  and 
theoretical  excess  volume  models  are  illustrated  in  Figure  5.5.  Because  there  was  no  data 
obtained  for  conversions  below  0.5,  it  is  anticipated  that  the  theoretical  models  will  be 
inaccurate  over  the  entire  conversion  range.  It  is  also  impossible  to  determine  if  the 
excess  volume  is  symmetric  over  the  entire  conversion  range.  It  can  be  seen  in  Figure  5.5 
that  neither  of  the  two  models  accurately  fit  the  experimental  data. 
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Figure  5.5.  Experimental  and  theoretical  excess  volume  versus  conversion 
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5.6.2  Conversion  Calculations 

The  one-  and  two-parameter  models  are  used  to  improve  the  conversion 
calculations  shown  in  Figure  5.6.  The  one-parameter  model  does  a better  job  of  matching 
the  gravimetric  conversion  than  the  two-parameter  model  except  between  0.8  and  0.9. 
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Figure  5.6.  Conversion  versus  time  for  different  models 
This  conclusion  could  change,  however,  with  the  addition  of  both  densitometer  and 

gravimetric  conversions  between  0 and  0.5.  The  excess  volume  may  turn  out  to  be 

symmetric  throughout  the  conversion  range,  allowing  the  two-parameter  model  to  be 

more  accurate.  More  experimentation  is  required. 

5.7  Polymerization  of  Butyl-Methacrylate 

5.7.1  Development  of  Models 

The  excess  volume  model  is  determined  by  taking  the  difference  between  the 
gravimetric  and  densitometer  conversion  at  different  conversions.  The  experimental  and 
theoretical  excess  volume  models  are  illustrated  in  Figure  5.7.  The  shape  of  the 
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experimental  excess  volume  is  not  symmetric  and  the  shape  does  not  match  any  of  the 
available  shapes  from  the  one-  and  two-parameter  models.  Both  the  one-  and  two- 
parameter  models  do  not  agree  with  the  experimental  data.  It  is  possible  that  the  models 
may  actually  decrease  the  accuracy  of  the  conversion  calculations. 
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Figure  5.7.  Experimental  and  theoretical  excess  volume 

5.7.2  Conversion  Calculations 

In  Figure  5.8,  the  gravimetric,  densitometer,  and  excess  volume  conversion 
calculations  are  plotted  versus  time.  In  this  example,  the  densitometer  conversion  is 
actually  more  accurate  than  the  conversions  determined  with  the  one-  and  two-parameter 
models.  Both  the  one-  and  two-parameter  models  diverge  from  the  gravimetric 
conversion  at  some  point  during  the  reaction. 
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Figure  5.8.  Conversion  calculation  for  different  models  and  methods 

5.8  Conclusions 

In  this  chapter,  it  has  been  shown  that  it  is  possible  to  determine  the  conversion  of 
an  emulsion  polymerization  reaction  by  measuring  the  density  of  the  emulsion  over  time. 
The  density  is  inversely  proportional  to  the  specific  volume  of  the  emulsion  reaction.  The 
conversion  of  the  reaction  can  then  be  measured  on-line  with  the  use  of  a densitometer. 

There  are  differences,  however,  between  the  conversion  determined  by 
gravimetric  analysis  and  the  conversion  calculated  by  the  densitometer.  This  difference 
can  be  modeled  through  the  use  of  an  excess  volume  model.  The  model  can  take  many 
mathematical  forms;  the  two  that  are  used  here  involved  both  a one-  and  two-parameter 
model. 

Deciding  which  model  to  use  is  the  critical  design  problem.  The  model  chosen 
must  be  able  to  guarantee  that  a conversion  root  will  be  found  inside  0 and  1 .0.  Different 
polymers  may  yield  different  results  and  may  require  that  a certain  model  is  used  in  order 


to  ensure  a correct  conversion  calculation. 


84 


It  has  been  shown  above  that  this  methodology  can  be  applied  to  several 
polymerization  reactions.  The  two-parameter  model  was  best  suited  to  styrene  while  the 
one-parameter  model  was  used  in  the  copolymerization  reaction.  It  was  even  illustrated 
with  butyl-methacrylate  that  sometimes  n ot  u sing  a n e xcess  v olume  m odel  i s t he  o nly 
practical  solution.  The  methodology  described  works  and  is  effective  for  an  emulsion 
polymerization  reaction. 


CHAPTER  6 

DESIGN  AND  CONTROL  OF  A CONTINUOUS  SAMPLING  AND  DILUTION 

SYSTEM 

6.1.  Introduction 

6.1.1  Motivation  Behind  the  CSDS 

In  recent  years,  real  time  monitoring  of  chemical  processes  has  become  very 
important  in  both  academia  and  industry.  On-line  monitoring  of  process  conditions  such 
as  concentration,  pressure,  or  temperature  is  highly  desirable  in  order  to  design  an 
effective  control  scheme  (Sacoto  1999).  On-line  monitoring  allows  for  improved  product 
quality  that  would  not  be  possible  through  off-line  techniques.  For  example,  real-time 
monitoring  allows  improvements  in  the  resolution  of  data. 

Real-time  measurements  are  practically  impossible  as  there  is  always  an 
associated  time  lag  with  the  steps  between  sampling  and  obtaining  the  desired  data. 
Regardless,  real-time  conditions  can  be  approached  if  the  correct  sensing  strategy  is 
utilized. 

There  are  four  different  types  of  strategies  that  were  considered  for  this  process: 
off-line,  on-line,  in-line,  andin-situ.  E ach  strategy  has  advantages  and  disadvantages. 
Off-line  requires  taking  a sample  and  performing  the  analysis  after  the  experiment  has 
concluded.  On-line  sampling  takes  a sample  and  analyzes  it  while  the  experiment  is  still 
running  while  in-line  sampling  introduces  a sensor  within  the  process  itself.  In-situ 
monitoring  is  the  least  invasive  of  the  four  techniques  as  it  places  a sensor  around  the 
experiment. 
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Although  on-line  sampling  methods  do  not  approach  real-time  measurements  as 
much  as  in-line  or  in-situ  measurements,  they  do  have  one  main  advantage.  On-line 
methods  can  use  existing  sensors  with  the  fewest  modifications  while  in-line  and  in-situ 
measurements  require  expensive  and  very  specialized  probes.  In  addition,  on-line  sensing 
is  very  similar  to  the  off-line  strategy  making  the  data  interpretation  similar.  Because  of 
the  reasons  stated  above,  the  on-line  strategy  was  employed  in  this  project. 

With  appropriate  s ampling  and  d ilution,  m uch  oft  oday’s  m ost  c ommonly  used 
sensing  equipment  can  be  modified  to  monitor  chemical  processes  on-line  (Sacoto  1999). 
In  the  area  of  concentrated  dispersions,  defined  as  highly  concentrated  solutions,  many 
researchers  have  attempted  to  monitor  composition  and  particle  properties.  For  example, 
emulsion  polymerization  researchers  have  often  measured  polymer  conversion  on-line 
through  the  use  of  a densitometer  (Canegallo,  et.  al,  1993;  Moritz,  et.  al,  1985)  or  free 
surfactant  concentration  from  surface  tension  (Schork  and  Ray,  1983).  Additionally  there 
is  a great  deal  of  research  in  the  area  of  on-line  measurements  of  particle  size  and  particle 
size  distribution  (Elicabe  and  Garcia-Rubio,  1988;  Kourti  et.  al.,  1984).  Many  of  these 
techniques  require,  or  would  be  greatly  improved  by,  some  type  of  dilution  system. 

This  chapter  reports  on  the  design,  implementation,  and  control  of  a continuous 
sampling  and  dilution  system  (CSDS)  developed  at  the  University  of  South  Florida  for 
on-line  monitoring  of  concentrated  dispersions. 

6.1.2  Issues  of  CSDS  Monitoring 

There  are  a number  of  challenges  associated  with  on-line  dilution  systems. 
Included  among  these  challenges  are:  representative  sampling,  time  delay  associated  with 
dilution,  and  sample  integrity.  This  section  covers  each  of  these  problems  in  turn. 
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Representative  sampling  involves  finding  a suitable  sampling  location,  the 
optimal  number  of  samples,  and  discovering  a suitable  amount  of  sample  required  to 
perform  the  analysis.  Locating  a suitable  sampling  location  is  crucial  since  a well-mixed 
vessel  is  not  necessarily  uniformly  distributed.  For  example,  in  a continuous  stirred  tank 
reactor  non-ideal  mixing  yields  dead  zones  while  in  concentrated  dispersions  the  size 
distribution,  surface  charge,  or  density  may  introduce  dead  zones  (Sacoto,  1999).  If  the 
sampling  location  is  chosen  in  a dead  zone  the  sample  is  not  representative  of  the  entire 
reactor.  One  solution  to  this  problem  is  to  sample  from  several  locations  but  the  added 
cost  and  complexity  associated  with  this  development  may  be  too  high. 

Depending  on  the  dilution  required,  on-line  sampling  may  require  large  amounts 
of  sample  which  would  adversely  affect  the  system,  i.e.  a process  disturbance.  If  the 
dilution  is  to  be  continuous,  then  large  amounts  of  waste  will  be  generated.  Because 
many  techniques,  light  scattering  or  chromatography  for  example,  require  large  dilutions 
this  becomes  a serious  problem.  For  UV-VIS  spectrometry  a dilution  factor  of  five  to  six 
orders  of  magnitude  is  required  (Sacoto,  1996).  Thus,  for  a given  dilution  step  10  to  100 
liters  of  diluent  is  required  for  every  milliliter  of  s ample  ( Sacoto,  1 999).  T he  diluted 
sample  does  not  return  to  the  experimental  process,  resulting  in  10  to  100  liters  of  waste 
per  milliliter  of  sample  analyzed. 

Time  delay,  or  transport  delay,  is  also  of  critical  importance.  If  the  dilution  is 
performed  in  steps  (as  in  the  CSDS)  then  each  dilution  step  increases  the  time  delay  and 
biases  the  interpretation  of  the  data.  The  time  delay  is  often  seen  to  increase  the  error  in 
the  estimation  of  particle  concentration  and  particle  size  distribution  (Elicabe  and  Garcia- 
Rubio,  1988).  There  are  two  major  effects  observed  as  the  time  delay  increases:  (1)  a 
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shift  in  time  relative  to  the  sample  populations  present  and  (2)  an  averaging  of  the  data 
over  the  measurement  interval  and  a loss  in  resolution. 

The  sample  integrity  must  be  preserved  from  the  time  the  sample  is  taken  from  the 
reactor  to  the  time  it  is  analyzed.  This  may  be  of  critical  importance  depending  on  the 
composition  of  the  dispersion.  For  example,  a change  in  the  continuous  phase  will  affect 
the  equilibrium  of  the  system;  an  alteration  in  the  chemical  potential  will  force  other 
phases  to  be  altered  to  keep  equilibrium.  Thus,  the  components  within  the  system  would 
be  redistributed  and  the  sample  analyzed  by  the  sensor  would  be  significantly  different 
from  what  originally  left  the  reactor.  For  example,  in  emulsion  polymerization  polymer 
particles  and  monomer  droplets  are  present,  but  when  diluted  the  monomer  droplets 
disappear.  The  particle  size  distribution  and  distribution  of  monomer  within  the  particles 
would  then  be  biased  and  the  measurement  inaccurate. 

6.1.3  Novelty  of  the  CSDS 

There  have  been  many  types  of  dilution  systems  developed  throughout  the  years. 
A system  developed  by  Nicoli  and  E lings  ( 1989)  consists  o f a s eries  o f d ilution  steps 
through  stirred  tanks  and  differential  pressure  for  sampling  as  shown  in  Figure  6.1.  A 
sample  is  taken  from  the  experimental  reactor  and  passed  into  the  first  CSTR  where  it  is 
diluted.  A second  sample  is  taken  from  the  first  CSTR  and  passed  into  a second  CSTR 
where  it  is  diluted  again.  The  process  is  repeated  n times  until  the  final  sample  is  sent  to 
the  sensor.  The  design  i s c ontinuous  i n n ature,  b ut  v ery  1 arge  a mounts  o f d iluent  a re 
required  and  there  are  time  distribution  problems  involved  with  t he  c ontinuous  s tirred 


tank  reactors  (CSTR). 
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Sample 


Figure  6.1.  Di 


n steps 


->  Sensor 


ution  system  proposed  by  Nicoli  et.  al.  that  can  be  visualized 
as  a sequence  of  CSTR  in  series. 

Ponnuswamy  and  Shah  (1986)  compiled  various  analytical  measurements  carried  both 
off-  and  on-line  to  monitor  polymer  quality  in  a batch  polymerization  reactor.  White 
(1994)  designed  a dilution  system  that  is  on-line  but  does  not  have  continuous  sampling 
and  dilution  using  various  solvent  and  sample  vessels  for  multiple  dilutions.  Most 
previous  work  on  dilution  systems  does  not  have  continuous  sampling  and/or  has  time 
distribution  problems  associated  with  it. 

6.2  Design  and  Implementation  of  CSDS 


6.2.1  General  design 

A continuous  sampling  and  dilution  system  developed  at  the  University  of  South 
Florida  consists  of  a tubing  network  of  interconnected  sampling  and  dilution  lines  that 
run  parallel  to  each  other  as  shown  in  Figure  6.2.  The  process  begins  when  a sample  is 
taken  from  the  reactor  and  mixed  with  a line  carrying  diluent.  The  sample  and  diluent  are 
then  allowed  to  mix  over  a designated  tube  length  while  static  mixers  are  added  to 
increase  radial  mixing  and  ensure  a homogeneous  bulk  concentration.  A sample  is  then 
taken  from  the  first  dilution  line  and  pumped  to  a second  dilution  channel  creating  a 
second  dilution  step.  The  dilution  process  is  repeated  n times. 
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Figure  6.2.  Schematic  of  the  CSDS  with  three  dilution  steps 
The  premise  behind  the  continuous  sampling  and  dilution  system  is  to  choose  the 
diluent  flow  rates,  Qj (t) , to  be  significantly  larger  than  the  sample  flow  rates,  q^\(t) , 
where  i is  the  number  of  dilution  steps.  Sacoto  (1999)  derived  the  following  steady  state 
equation  for  concentration 

c (t)  = Q-i(Qgj-i(0 
QiiO  + q,- 1(0 

where  Q (t)  corresponds  to  the  overall  sample  concentration  at  the  ib  step.  It  should  be 
noted  that  if  Q,  ( t ) > qt  (t)  it  will  lead  to  very  high  dilution  ratios.  However,  this 
equation  only  analyzes  one  dilution  step;  for  multiple  steps  these  equations  need  to  be 
combined  to  determine  the  end  concentration. 

There  are  many  advantages  to  using  the  continuous  sampling  and  dilution  system 
as  opposed  to  other  dilution  systems  such  as  the  system  developed  by  Nicoli  and  Elings 
(1989).  The  CSDS,  unlike  some  previously  designed  systems,  does  not  suffer  from  age 
distribution  problems  because  the  CSDS  can  be  modeled  as  several  plug  flow  reactors  in 


series. 
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The  CSDS  is  a modular  system;  extra  dilution  steps  may  be  added  or  removed 
easily.  If  the  dispersion  is  highly  concentrated,  it  may  be  necessary  to  add  a step.  This 
has  the  adverse  effect  of  increasing  the  dilution  time,  i.e.,  the  time  from  when  the  sample 
is  taken  to  when  the  sample  is  analyzed  by  the  spectrometer.  However,  even  with  several 
dilution  steps,  the  time  lag  is  on  the  order  of  5 minutes. 

Next,  the  dilution  times  are  short,  roughly  5 or  6 minutes  (Sacoto,  1999),  which 
allows  faster,  more  accurate  control.  Finally,  the  CSDS  can  be  used  continuously,  which 
makes  it  a potentially  powerful  tool  for  process  control. 

6.2.2  The  Dynamic  Model  of  the  CSDS 

In  order  to  implement  advanced  control  strategies  (such  as  a Smith  predictor  or 
model  predictive  control)  an  accurate  model  of  the  system  must  be  developed.  A first 
principles  model  is  derived  below  and  is  used  for  simulation  and  controller  design.  The 
derivation  begins  by  performing  a mass  balance  on  the  first  stage  of  the  dilution  system 

= &(t)C  + qQ(t)C0(t)-qi  (t)Cx  (t)  - 0lo  (0  Q (0  (6-1) 

at 

There  is  no  product  in  the  dilution  line,  thus  C = 0.  Q\  is  defined  as 

0,o(  O = 0i(O  + 0o(O-tf,(O 


Plugging  this  expression  into  equation  (6.1)  yields 

= ?0«)Q(o  - «,  to c, (o  - (am +?„(<)  - «,m)c, m 


The  terms  ql(t)Cl(t)  cancel  each  other  out  leaving  the  following  equation 


d{VxCx) 


9o(OQ(o-(a(o+9o(o)cI(o 


dt 
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Because  the  volume  is  constant  in  the  dilution  system,  it  can  be  factored  out  of  the 
integral  and  simplified  to  result  in 


jgS>,&Wc;(o-(Q(0+g  °(0)cl(o 

dt  V. , 0 F,  1 


(6.2) 


Equation  (6.2)  is  then  checked  against  the  steady  state  model  developed  by  Sacoto  (1999) 
shown  below 

rfA._  Q(Qgo(0 

,U  0(O+*o(O 


The  derivative  in  equation  (6.2)  is  set  equal  to  zero  in  order  to  calculate  a steady  state 
solution  and  rearranged  to  give 


go(0 

K 


C0(  0 = 


(gi(O  + g0(O) 


cx(t) 


(6.3) 


The  volumes  cancel  each  other  and  equation  (6.3)  can  be  rearranged  to  yield 

Q(0?o(0 
lU  Q(0+9o(0 

which  proves  that  the  dynamic  equation  derived  is  equivalent  to  the  steady  state  equation 
as  reported  by  Sacoto  (1999).  A generic  dilution  step  model  is  shown  below  and  is 
similar  in  structure  to  equation  (6.2) 

(0-(Q(Q+gM(Q)C(0  (64) 

dt  Vf  V; 

where  i is  the  dilution  step. 

6.2.3  Implementation 

The  CSDS  prototype  has  been  implemented  as  depicted  in  Figure  6.2.  There  are 
four  main  parts  to  the  system:  the  sampling  system,  the  pumping  system,  the  tubing 
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network,  and  the  sensor.  Each  of  these  sections  may  be  modified  depending  on  the 
application  though  for  this  discussion  the  most  applicable  scheme  developed  at  the 
University  of  South  Florida  is  implemented. 

The  reactor  sampling  is  carried  out  by  a sample  recycle  stream  from  which  a 
sample  is  taken  and  sent  to  the  CSDS.  The  most  important  design  features  here  are  the 
tubing  selection  and  pump  speed  which  must  be  chosen  to  minimize  the  sample  and 
prevent  clogging. 

There  are  two  different  pumps  used  to  operate  the  dilution  system:  sample  pumps 
and  dilution  pumps.  Both  pumps  use  the  same  model  pump  heads  and  drivers;  peristaltic 
quick  load  made  by  Cole-Parmer  (Model:  7021-24)  and  a modular  drive  by  Cole-Parmer 
(Model  7553-70)  respectively.  The  advantage  to  using  two  separate  pumps  is  that  the 
dilution  ratio  can  be  varied  greatly  because  the  sample  and  dilution  lines  can  be 
manipulated  independently  of  one  another. 

The  flow  rates  of  the  sample  and  dilution  lines  depend  on  the  RPMs  of  the 
respective  pumps  and  the  inner  diameter  of  the  tube.  Therefore  tube  diameter  is  a critical 
design  criterion  and  the  pump’s  RPMs  may  be  manipulated  to  maintain  a desired 
concentration  of  sample  in  the  system.  The  greater  the  difference  between  the  inner 
diameter  of  the  sample  line  and  the  inner  diameter  (ID)  of  the  dilution  line,  the  larger  the 
dilution  factor  (assuming  constant  RPM  for  both  the  sample  and  dilution  lines). 
However,  the  larger  the  ID  of  the  dilution  lines,  the  more  diluent  required.  Conversely, 
the  smaller  the  sample  line  ID,  the  greater  the  chance  that  the  line  will  clog.  An  optimal 
design  must  be  developed  for  the  particular  application  that  takes  these  two  factors  into 
account.  Some  of  the  other  important  design  features  of  the  CSDS  tubing  system  that  are 
desired  include:  minimum  residence  time  in  the  CSDS,  minimum  sample  volume  taken 


94 


from  reactor,  avoidance  of  clogging  (caused  by  high  RPMs  and  small  tube  diameters), 
and  prevention  of  temperature  gradients. 

The  tubing  material  must  also  be  considered  as  the  interaction  of  the  tubing 
material  with  the  solvent  may  have  possibility  to  absorb  part  of  the  sample  on  the  tubing. 
For  example,  in  polystyrene  latex  reactions,  styrene  adsorbs  on  the  walls  of  many  tubing 
materials  (Sacoto,  1999).  Suitable  materials  must  be  found  depending  on  the  solvent  and 
substance  analyzed. 

A UV-VIS  spectrometer  built  by  Ocean  Optics  was  chosen  for  this  project.  The 
spectrometer  measures  the  absorbance  spectra  from  190nm  to  820nm  and  was  connected 
to  a desktop  PC  for  analysis.  A flow  through  cell  was  used  to  interface  the  spectrometer 
with  the  C SDS.  T he  C SDS  final  dilution  stream  is  connected  to  the  this  cell  and  the 
spectroscopic  analysis  is  performed  continuously. 

6.3  Automation  and  Control  of  the  CSDS 

The  main  purpose  of  the  CSDS  is  to  dilute  the  product  so  that  the  absorbance,  i.e. 
concentration,  is  within  a specified  window  so  spectroscopic  analysis  can  be  performed 
on  the  sample.  In  this  system,  the  sample  pump’s  rpm  is  kept  at  a constant  rate;  thus 
changing  the  dilution  pump’s  rpm  will  alter  the  dilution  ratio,  and  likewise  the 
absorbance.  The  automation  of  the  CSDS  consists  of  getting  the  computer,  pumps,  and 
spectrometer  communicating  with  one  another  in  order  to  produce  the  correct  absorbance. 

LabVIEW  and  National  Instruments  data  acquisition  boards  are  used  to  interface 
the  dilution  system  with  the  computer.  Along  this  line  several  Vis  (virtual  instruments) 
had  to  be  written  and  developed  in  LabVIEW  to  communicate  with  the  pumps  and 
spectrometer  and  the  overall  scheme  and  flow  of  information  are  illustrated  in  Figure  6.3. 
There  were  four  Vis  that  had  to  be  developed:  the  pump,  dilution  controller,  Ocean 
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Optics  (spectrometer),  and  post-processing  Vis.  The  process  begins  with  the  user 
specifying  an  absorbance  value.  The  dilution  controller  calculates  a dilution  pump 
voltage  based  on  the  equation  for  a PI  controller.  This  voltage  setpoint  is  passed  to  the 
pumps  via  the  pump  VI  which  p asses  t he  v oltage  t hrough  t he  c orrect  c hannels  o n t he 
DAQ  (data  acquisition)  board  to  the  pumps.  The  dilution  pump’s  RPMs  change  as  a 
function  of  the  new  voltage  passed  to  them,  thus  the  concentration  (and  absorbance)  of 
product  changes  as  shown  in  equation  (6.4).  The  diluted  sample  is  then  passed  into  the 
flow  through  cell  where  the  spectral  data  is  then  read  from  the  spectrometer  by  the  Ocean 
Optics  VI  and  sent  to  the  Post  Processing  VI.  Here  the  spectral  data  and  dilution  ratio  are 
used  to  calculate  the  product  concentration  and  absorbance  value.  The  absorbance  is  sent 
back  to  the  dilution  controller  to  make  any  further  adjustments. 


iolyner 

Coaceatraitioi 

Figure  6.3.  Schematic  of  control  scheme  in  Lab  VIEW 
A PI  controller  was  coded  into  LabVIEW  and  is  used  for  the  dilution  controller. 

The  proportional  and  integral  parameters  were  determined  by  doing  a step  response  and 

using  the  ITAE  settings  and  will  be  explained  in  greater  detail  later.  In  order  to  help  the 

controller,  anti-reset  windup  was  added  to  the  design.  There  are  two  parts  to  this: 


96 


installing  a maximum  upper  limit  on  the  error  signal  sent  to  the  controller  and  allowing 
the  user  to  clear  out  the  integral  error  present  as  shown  by  Mattem  (1993). 

The  user  interface  is  shown  in  Figure  6.4.  The  user  has  the  ability  to  specify  the 
absorbance  setpoint,  change  the  voltage  (RPM)  of  the  sample  pumps,  change  the 
controller  parameters  (proportional,  integral,  and  bias),  and  specify  the  sample  time.  The 
user  also  has  control  over  the  anti-reset  windup  by  changing  the  upper  and  lower  error 
limits.  The  PI  controller  must  first  be  initialized  and  then  placed  in  run  mode  for  the 
system  to  function  properly.  Finally,  the  data  can  be  saved  to  a spreadsheet  if  so  desired. 
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Figure  6.4.  User  interface  for  controlling  the  CSDS  in  LabVIEW 
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6.4  Experimental  and  Simulated  Responses  of  the  CSDS 
6.4.1  CSDS  Open-Loop  Response 

The  CSDS  was  run  in  open-loop  mode  (i.e.  the  controller  was  turned  off)  and  a 
step  change  made  in  the  dilution  pump’s  voltage  from  5V  to  7V.  The  absorbance  value 
was  recorded  as  a function  of  time  and  the  result  is  illustrated  in  Figure  6.5.  There  is  a 
short  time  delay  ( 9 ) of  25s  followed  by  a sharp  decrease  in  absorbance.  The  increase  in 
voltage  to  the  dilution  pump  causes  the  RPMs  to  increase  and  likewise  the  diluent  flow 
rate  increases.  This  rise  in  diluent  flow  rate  produces  a decreased  concentration  of  final 
sample  analyzed  in  the  spectrometer  and  thus  a lower  absorbance  value. 


Unfiltered 

-Filtered 


Figure  6.5.  Open-loop  response  of  the  CSDS  to  a step  change  from  5V  to 
7V  in  the  dilution  pump’s  voltage. 

6.4.2  CSDS  Closed-Loop  Response 

Figure  6.5  can  be  approximated  as  a first  order  plus  time  delay  (FODT)  system 
and  the  parameters  K,  9,  and  r can  be  easily  determined  from  such  experimental  response 
data  (Seborg  et.  ai,  1989)  where  K is  the  gain,  9 is  the  time  delay,  and  r is  the  time 
constant.  One  approach  to  developing  controller  relationships  is  based  on  a performance 
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index  that  takes  into  account  the  entire  closed-loop  system.  T he  i ntegral  o f t he  t ime- 
weighted  absolute  error  (ITAE)  is  one  method  based  on  the  equation 

ITAE  = J®  t\e(t)\dt 

where  PI  controller  correlations  have  been  developed  that  minimize  this  error.  The  ITAE 
formula  penalizes  e rrors  t hat  e ndure  f or  1 ong  p eriods  o f t ime  a nd  p referred  o ver  o ther 
methods  because  it  results  in  conservative  controller  settings.  For  this  project,  the  type  of 
input  considered  was  a set  point  change  which  results  in  the  following  equations  for  a PI 
controller 

KKC  = 0.586(0/ r)-0'916 
T/Ti  =1.03  -0.1 65(0/ r) 

where  Kc  is  the  controller  gain  and  r,  is  the  integral  time.  Using  the  step  response  in 
Figure  6.5  the  following  values  for  controller  parameters  were  calculated 

Kc  - -3.55  r,  = 33.65 

It  i s w orth  n oting  t hat  t he  g ain  0 f t he  c ontroller  i s n egative  b ecause  t his  is  an  inverse 
response  system  meaning  that  an  increase  in  the  diluent  flow  rate  will  decrease  the 
absorbance  (concentration)  of  the  final  sample. 

Phenol  red  solution  combined  with  a buffer  solution  (pH  7.4)  was  used  to  test  the 
dilution  s ystem.  The  sample  pump’s  rpms  were  kept  at  a constant  level  (100  rpm)  in 
order  to  create  a constant  sample  environment.  The  ITAE  tuning  parameters  were  entered 
into  the  system  via  the  user  interface  in  Lab  VIEW.  The  controller  was  initialized  and  a 
step  change  was  made  in  the  absorbance  from  0.7AU  to  0.9AU.  In  Figure  6.6  the  closed- 
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loop  response  of  the  CSDS  is  illustrated.  The  system  reaches  the  new  steady  state  set 
point  in  less  than  4 minutes  and  produces  a stable,  smooth  response. 
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Figure  6.6.  Closed-loop  response  of  absorbance  with  ITAE  tunings 
parameters 

The  response  of  the  input  (dilution  voltage)  to  the  CSDS  was  also  recorded  versus  time 
and  is  shown  in  Figure  6.7.  The  voltage  begins  at  roughly  7.5V  and  decreases  sharply  to 
a value  around  6V.  The  voltage  decreases  because  the  set  point  of  the  absorbance 
increases  and  the  CSDS  is  an  inverse  response  system. 


time  [s] 


♦ Dilution  voltage 


Figure  6.7.  Closed-loop  response  of  dilution  voltage  (input)  with  ITAE 
tuning  parameters 
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6.4.3  CSDS  Simulation 

The  dynamic  model  of  the  dilution  system  was  coded  into  Matlab  and  Simulink 
and  a S-function  written  in  Matlab  entitled  csds_control.m  and  is  shown  in  Figure  6.8. 
Csds_control  contains  the  differential  equations  from  the  dynamic  model  that  are 
responsible  for  calculating  the  changes  in  c oncentration  ( absorbance)  f or  e ach  d ilution 
step.  Transport  delays  are  also  included  in  the  model,  but  are  currently  static  delays 
despite  the  fact  that  the  RPMS  of  the  pump  will  vary  causing  changing  time  delays. 


Figure  6.8.  Diagram  of  the  Simulink  model  used  for  simulations 
A model  for  the  relationship  between  the  pump  speed  and  time  delay  must  be 

established  in  order  to  be  used  in  MATLAB  and  Simulink.  A simple  linear  model  was 

proposed  to  model  the  pumps  and  the  parameters  estimated  using  graphs  similar  to  that 

illustrated  in  Figure  6.9.  Each  line,  dilution  and  sample  has  its  own  corresponding 

equation  that  is  used  to  model  the  pump. 
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Figure  6.9.  Example  of  pump  calibration  curve  used  for  CSDS  simulation 
The  models  for  the  dilution  and  sample  pumps  were  obtained  and  coded  into 

MATLAB.  A simulation  of  t he  C SDS  w as  p erformed  u sing  t he  s ame  e xact  s et  p oint 

change  and  controller  parameters  as  discussed  in  section  6.4.2.  As  seen  in  Figure  6.10, 

the  simulated  response  agrees  reasonably  well  with  the  experimental  data.  Both  of  the 

responses  exhibit  the  same  time  delay  and  have  similar  rise  times. 


Figure  6.10.  Theoretical  and  experimental  responses  of  first  dilution  step  to 
setpoint  change 
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The  oscillatory  nature  of  the  two  responses  exhibit  the  same  trends,  yet  have  distinct 
differences.  For  example  the  simulated  response  leads  the  experimental  response  both  in 
rise  time  and  in  the  changes  in  oscillation.  This  difference  could  be  due  in  part  to  the 
static  time  delay  used  in  the  simulation  while  the  experimental  system  has  a dynamic  time 
delay. 

Further  agreement  can  be  seen  in  Figure  6.11  that  compares  the  theoretical  and 
experimental  dilution  pump  responses.  The  shapes  of  the  two  curves  are  almost  identical 
being  only  slightly  offset  from  one  another.  The  encouraging  aspect  is  that  the  slopes  are 
the  same  for  the  entire  time  frame  involved. 


Figure  6.1 1.  Theoretical  and  experimental  dilution  voltage  responses 

6.5  Conclusions  and  Future  Work 

There  are  several  places  in  the  dilution  system  where  more  work  needs  to  be  done. 
First,  a flowmeter  must  be  installed  and  verified  which  will  provide  more  accurate  models 
for  the  pumps  in  terms  of  both  transient  and  steady-state  behavior.  These  more  accurate 
models  will  aid  both  experimental  and  theoretical  work  by  improving  controller  tuning 
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parameters.  The  most  important  aspect  of  the  flowmeter,  however,  will  be  the  added 
ability  to  measure  the  dilution  ratio  instead  of  using  models  to  calculate  it.  This  will  lead 
to  a more  accurate  determination  of  the  concentration. 

As  more  dilution  steps  are  added,  the  time  delay  will  increase.  A more  advanced 
control  scheme,  such  as  a Smith  predictor  or  predictive  control,  may  be  necessary.  It  is 
believed  that  a Smith  predictor  or  predictive  control  scheme  would  greatly  enhance  the 
ability  of  the  controller  to  react  to  both  disturbances  and  changes  in  set  point. 


CHAPTER  7 

CONCLUSIONS  AND  FUTURE  WORK 

7.1  Conclusions 

This  dissertation  presents  an  in-depth  exploration  of  theoretical  robust  control  as 
well  as  the  design  and  implementation  of  subsystems  for  an  emulsion  polymerization 
reactor. 

Chapters  2 through  4 investigate  theoretical  process  control  challenges, 
specifically  theoretical  robust  control  and  disturbance  predictive  control.  Chapter  2 
presents  a revised  scheme  of  the  critical  direction  theory  that  is  now  applicable  to 
nonlinear,  non-convex  systems.  In  addition,  analytical  expressions  are  determined  for  the 
minimizing  amplitude  for  a variety  of  nonlinear  systems.  Chapter  3 analyzes  the 
robustness  of  predictive  controllers  using  the  parametric  stability  margin.  In  particular,  a 
systematic  algorithm  is  developed  that  finds  predictive  controller  tuning  parameters  that 
make  the  controller  robustly  stable  for  a given  uncertainty.  Chapter  4 develops  a 
disturbance  predictive  control  law  that  is  applicable  to  systems  that  have  large 
disturbances  present  such  as  a wastewater  treatment  plant.  It  was  shown  that  for  linear 
plants  this  control  law  does  an  excellent  job  of  disturbance  rejection  but  performs  poorly 
for  highly  nonlinear  plants. 

Chapter  5 discusses  the  addition  of  an  excess  volume  model  to  conversion 
calculations.  The  excess  volume  model  increases  the  numerical  complexity  of  the 
analysis  but  greatly  improves  the  conversion  calculations  for  a variety  of  systems. 
However,  the  numerical  issues  were  analyzed  and  a solution  obtained.  Finally,  Chapter  6 
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extends  the  work  on  the  CSDS  that  was  begun  at  the  University  of  South  Florida.  The 
CSDS  was  automated  and  controlled  via  LabVIEW  and  a dynamic  model  was  developed 
for  the  system.  The  control  scheme  implemented  keeps  the  absorbance  value  of  the 
sample  at  a specified  threshold  required  for  spectroscopic  measurements. 

7.2  Future  Work 

The  nonlinear  Nyquist  robust  stability  margin  developed  in  Chapter  2 needs  to  be 
applied  to  relay  based  PI  controllers  and  perhaps  extended  to  include  controller  synthesis. 
Another  a rea  o f r esearch  i nterest  i s t he  s tability  a nalysis  o f i nfinite  h orizon  p redictive 
controllers  developed  by  Cannon  and  Kouvaritakis  (2000).  Work  is  currently  proceeding 
on  developing  the  algorithms  involved  and  determining  if  the  infinite  horizon  increases  or 
decreases  the  robustness  of  the  controller.  The  disturbance  predictive  controller 
developed  in  Chapter  3 must  be  modified  so  it  is  effective  for  nonlinear  systems.  This 
could  be  done  through  an  adaptive  control  scheme  or  more  effectively  through  the  use  of 
a nonlinear  disturbance  predictive  controller. 

Finally  the  subsystems  developed  in  Chapters  5 and  6 need  to  be  implemented  on 
an  emulsion  polymerization  reactor.  In  conjunction  with  a multi-angle,  multi-wavelength 
spectrometer  the  particle  size  distribution  can  be  obtained  on-line.  Significant  advances 
in  the  control  of  the  particle  size  will  be  possible  once  the  test  bed  is  implemented. 


APPENDIX  A 

PROOF  OF  EQUATION  (5.15) 

Proof.  The  following  proof  is  derived  using  root  analysis  for  a cubic  equation. 
The  general  form  of  any  cubic  equation  is 

3 2 

z +«2Z  +ci\z  + ao  = 0 

and  define  q and  r as 

q=\a'~\a*  (A-1) 


r = ^{ala2~3ao)-^a2 


(A. 2) 


It  can  be  shown  that  if 


D = q3  +r2  >0 


(A.3) 


there  is  one  real  root  and  two  complex  conjugate  roots.  The  cubic  equation  for 
conversion  using  the  Redlich-Kister  model  is 


x3  + 


^ 35  7 f A + v]nft0/„  vno/„  B^\  f Vq0/o  v(x ) 


2 B 


\x2  + 


2 B 


x + 


2 B 


= 0 


and  using  equations  (A.l),  (A.2),  and  (A.3)  the  discriminant  D becomes 


D = — 


J (65vioqo/0  - 65vqo/0  + 5 2 + A2  )3 

46656  B6 


+ 


(3B-A),. 

(^  + vioo% 

2452 


- v0o/o  - B)  + 


vp%  ~ v(*) 

45 


+ 


(35-  X)3 
21653 


(A.4) 


The  issue  in  determining  if  D > 0 involves  analyzing  the  cubic  term  in  (A.4)  to  see  if  it  is 
positive  or  negative.  Therefore,  define 
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6Bvm%  - 6Bv()%  + B2  + A2  <0 


and  rearrange  to  yield 


A2  +B2 
6B 


' < V0%  V100% 


The  above  equation  matches  equation  (5.15)  and  the  proof  is  concluded. 
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